This notebook takes the data generated by the script in /Ex_Quartic, which describes the newly developed method, and compares it existing classification algorithms. This multi-objective problem was selected as a simple illustration, as it only takes 2 inputs and provides 2 outputs, but has the complication of a local minimum and maximum besides the Pareto optimum.
The pattern recognition algorithms are: * Support Vector Machines: example supervised learning problem * Gaussian mixture models: example unsupervised learning problem
These are compared to the process of characterizing the optimal inputs using Gaussian Processes refined by adaptive sampling.
For ease of calculation, the iterative refinement of the GP models are not included here, nor are the calculations of the conditional probabilities. The SVM and GMM methods are tested with (1) the data after finding the Pareto front (before the new adaptive sampling process), (2) the data after adaptive sampling refinement, and (3) the data after the Pareto front search with additional random sampling to see if the new sampling method is useful for these other processes as well.
Comparison among the supervised methods is done by looking at the error rate across the entire space. Since the test function is a simple polynomial, the solution of what is acceptable can be found explicitly. For the GP method, since acceptance is defined by a probability, the error rate is weighted by the probability of acceptance. For the SVM method, the error rate is simply the number of incorrectly categorized points divided by the total number of test points.
In addition to the accuracy comparison, a comparison of the importance ranking metric developed for using GP in classification problems will be compared to Shapley values. As with the accuracy comparison, this will be tested against the data before refinement, after refinement, or after a random sample.
Clear the workspace and define the functions.
# Setup
rm(list = ls())
# Visualization
library(dplyr)
library(ggplot2)
library(patchwork)
# Parallel processing
library(parallel)
library(doParallel)
# Gaussian processes
library(GPareto)
library(DiceKriging)
library(DiceOptim)
# Optimization
library(GA)
# Gaussian Mixture Models
library(mclust)
__ ___________ __ _____________
/ |/ / ____/ / / / / / ___/_ __/
/ /|_/ / / / / / / / /\__ \ / /
/ / / / /___/ /___/ /_/ /___/ // /
/_/ /_/\____/_____/\____//____//_/ version 5.4.7
Type 'citation("mclust")' for citing this R package in publications.
# Support Vector Machines
library(e1071)
Relevant functions
The quartic functions to optimize.
# Set (x1, x2) on range of [0, 5]
# Objective functions are based on the quadratic functions used previously, but with an additional local minimum
f1 = function(x1, x2){
return(20*(x1 - 0.75)^2 + 190 + 11.58*x2^4 - 115.85*x2^3 + 383.13*x2^2 - 463.50*x2)
}
# The second objective function is also partially rotated so the local optimum is not perfectly aligned
f2 = function(x1, x2){
# Remap both variables: rotate 30 degrees counterclockwise
ang = -pi/24
n1 = x1*cos(ang) - x2*sin(ang)
n2 = x1*sin(ang) + x2*cos(ang)
# return((x1 - 2.5)^2 + 80 + 1.778*x2^4 - 20*x2^3 + 78.573*x2^2 - 124.664*x2)
return((n1 - 2.5)^2 + 80 + 1.778*n2^4 - 20*n2^3 + 78.573*n2^2 - 124.664*n2)
}
# Using GPareto: Need inputs and outputs as vectors/matrices not as dataframes. Coarse initial design
fun = function(x){
x1 = x[1]; x2 = x[2]
return(c(f1(x1, x2), f2(x1, x2)))
}
Normalization-related functions
# Normalized objective functions
n.obj = function(GPar.data, GPar.front){
# Given dataframes that describe the entire dataset and the front, find the normalized (x,y)
# Objective functions are named 'f1' and 'f2'
# Normalize the objective outputs so that the utopia point is (0,0) and the nadir point is (1,1)
f1.up = GPar.front$f1[which.min(GPar.front$f2)]
f2.up = GPar.front$f2[which.min(GPar.front$f1)]
GPar.data$f1.norm = (GPar.data$f1 - min(GPar.front$f1))/(f1.up - min(GPar.front$f1))
GPar.data$f2.norm = (GPar.data$f2 - min(GPar.front$f2))/(f2.up - min(GPar.front$f2))
return(GPar.data)
}
# Calculate the normalized distance
n.dist = function(f1.norm, f2.norm, GPar.front){
# Given the normalized coordinates (f1.norm, f2.norm) and the Pareto frontier estimate,
# find the distance along the constant f2/f1 ratio line
# Determine the two points on the Pareto front that define the relevant segment
GPar.front$theta = atan(GPar.front$f2.norm / GPar.front$f1.norm)
if(f1.norm < 0){f1.norm = 0}
if(f2.norm < 0){f2.norm = 0}
ratio = atan(f2.norm/f1.norm)
# Check if the angle is the same as a point on the Pareto front
if(any(abs(ratio - GPar.front$theta) < 1e-5)){
pos = which.min(abs(ratio - GPar.front$theta))
Par.x = GPar.front$f1.norm[pos]
Par.y = GPar.front$f2.norm[pos]
} else{ # Otherwise, two points are needed for linear interpolation
# Break the dataframe into theta above and below
Par.above = GPar.front[GPar.front$theta - ratio > 0,]
Par.below = GPar.front[GPar.front$theta - ratio < 0,]
# Find the point closest to the angle
pos.above = which.min(abs(ratio - Par.above$theta))
pos.below = which.min(abs(ratio - Par.below$theta))
# Linear interpolation
ln.x = c(Par.above$f1.norm[pos.above], Par.below$f1.norm[pos.below])
ln.y = c(Par.above$f2.norm[pos.above], Par.below$f2.norm[pos.below])
slp = diff(ln.y)/diff(ln.x)
# Find the point on the segment with the same angle, ie. the same ratio.
# Solving with this constraint has analytical solution:
Par.x = (ln.y[1] - slp*ln.x[1]) / (f2.norm/f1.norm - slp)
Par.y = slp*(Par.x - ln.x[1]) + ln.y[1]
}
# Linear distance to the front is the difference between distances to the origin
dist = sqrt(f1.norm^2 + f2.norm^2) - sqrt(Par.x^2 + Par.y^2)
return(dist)
}
Gaussian process parameter tuning
fill.sample.mod = function(GPar.data, input.name, output.name){
# Calculate the GP model to use.
# Using the km function, but applies checks on the system to make sure that
# the model uncertainty matches expectations based on GP, ie. it did not
# fail to converge.
# Based on testing, the model is bad when the 10% percentile and 90% percentile
# of the standard deviation are of the same order of magnitude. This is easiest
# checked if the difference between the 10th and 90th percentile
# is larger than the difference between the 25th and 75th.
pt10 = 1; pt90 = 1; pt25 = 1; pt75 = 1
while(log10(pt90/pt10) <= log10(pt75/pt25)){
mod.out = km(design = GPar.data[, input.name], response = GPar.data[, output.name],
covtyp = 'gauss', # Gaussian uncertainty
optim.method = 'gen', # Genetic algorithm optimization
control = list(trace = FALSE, # Turn off tracking to simplify output
pop.size = 50, # Increase robustness
max.generations = 400), # Some convergence issues
nugget = 1e-6, # Avoid eigenvalues of 0
)
# Randomly sample 200 points from the search space
pt = 200; i = 1
lims = range(GPar.data[,input.name[i]])
samp = data.frame(runif(n = pt, min = lims[1], max = lims[2]))
for(i in 2:length(input.name)){
lims = range(GPar.data[,input.name[i]])
samp[,i] = runif(n = pt, min = lims[1], max = lims[2])
}
names(samp) = input.name
# Find model output to find the percentile ranks for this iteration
res = predict(object = mod.out, newdata = samp, type = 'UK')
pt10 = quantile(res$sd, 0.10); pt90 = quantile(res$sd, 0.90)
pt25 = quantile(res$sd, 0.25); pt75 = quantile(res$sd, 0.75)
}
return(mod.out)
}
Comparison of the boundaries to show how changing the acceptance criteria changes the shape of the near-Pareto set. Showcases the robustness to different selection criteria, indicating flexibility in the design objectives.
## Loading
# Load datasets for obtaining the refined probability functions
data.delta = read.csv('../Ex_Quartic/GPar_Accept_Delta1.csv')
data.cutof = read.csv('../Ex_Quartic/GPar_Accept_Threshold.csv')
data.radan = read.csv('../Ex_Quartic/GPar_Accept_Radius.csv')
# Load datasets prior to refinement
data.paret = read.csv('../Ex_Quartic/GPar_all_start.csv')
data.paret$rad = sqrt(data.paret$f1.norm^2 + data.paret$f2.norm^2)
# Compare to the estimate of the Pareto frontier
GPar.front = read.csv(file = '../Ex_Quartic/GPar_fnt_start.csv')
# # Add a random samples to simulate the effect of sample size rather than adaptive sampling
# nsamp = max(nrow(data.delta), nrow(data.cutof), nrow(data.radan)) - nrow(data.paret)
# data.rando = data.frame(x1 = runif(n = nsamp, min = 0, max = 5),
# x2 = runif(n = nsamp, min = 0, max = 5))
# data.rando$f1 = f1(x1 = data.rando$x1, x2 = data.rando$x2)
# data.rando$f2 = f2(x1 = data.rando$x1, x2 = data.rando$x2)
# # Fill in the remaining calculations: normalized outputs, distance, theta, order
# data.rando = n.obj(GPar.data = data.rando, GPar.front = GPar.front)
# cl <- makeCluster(2)
# registerDoParallel(cl)
# dist = foreach(row = 1:nrow(data.rando)) %dopar%
# n.dist(f1.norm = data.rando$f1.norm[row], f2.norm = data.rando$f2.norm[row], GPar.front = GPar.front)
# stopCluster(cl)
# data.rando$dist = unlist(dist)
# data.rando$rad = sqrt(data.rando$f1.norm^2 + data.rando$f2.norm^2)
# data.rando$theta = atan(data.rando$f2.norm/data.rando$f1.norm)*180/pi*10/9
# data.rando$order = seq(from = max(data.paret$order) + 1, to = max(data.paret$order) + nsamp, by = 1)
# data.rando = rbind(data.paret[, names(data.paret) %in% names(data.rando)], data.rando)
#
# # Store the random data
# write.csv(data.rando, file = 'GPar_Random.csv')
data.rando = read.csv(file = 'GPar_Random.csv')
##
# Grid of the relevant region to visualize and compare
lower = c(0, 0); upper = c(5,5); grid.sz = 100
fine.grid = expand.grid(x1 = seq(from = lower[1], to = upper[1], length.out = grid.sz),
x2 = seq(from = lower[2], to = upper[2], length.out = grid.sz))
fine.grid = fine.grid[,1:2]
names(fine.grid) = c('x1', 'x2')
# Calculate the actual results
fine.grid$f1 = f1(x1 = fine.grid$x1, x2 = fine.grid$x2)
fine.grid$f2 = f2(x1 = fine.grid$x1, x2 = fine.grid$x2)
fine.grid = n.obj(GPar.data = fine.grid, GPar.front = GPar.front)
fine.grid$dist = NaN
for(row in 1:nrow(fine.grid)){
fine.grid$dist[row] = n.dist(f1.norm = fine.grid$f1.norm[row],
f2.norm = fine.grid$f2.norm[row],
GPar.front = GPar.front)
}
fine.grid$ang = atan(fine.grid$f2.norm/fine.grid$f1.norm)*180/pi*10/9
fine.grid$rad = sqrt(fine.grid$f1.norm^2 + fine.grid$f2.norm^2)
##
# Normalized distance
fine.grid$delta = 0
fine.grid$delta[fine.grid$dist <= 1] = 1
mod.dist.paret = fill.sample.mod(GPar.data = data.paret, input.name = c('x1', 'x2'), output.name = 'dist')
mod.dist.adapt = fill.sample.mod(GPar.data = data.delta, input.name = c('x1', 'x2'), output.name = 'dist')
mod.dist.rando = fill.sample.mod(GPar.data = data.rando, input.name = c('x1', 'x2'), output.name = 'dist')
res = predict(object = mod.dist.paret, newdata = fine.grid[,c('x1', 'x2')], type = "UK")
fine.grid$prob.delta.paret = pnorm(q = 0, mean = res$mean - 1, sd = res$sd)
res = predict(object = mod.dist.adapt, newdata = fine.grid[,c('x1', 'x2')], type = "UK")
fine.grid$prob.delta.adapt = pnorm(q = 0, mean = res$mean - 1, sd = res$sd)
res = predict(object = mod.dist.rando, newdata = fine.grid[,c('x1', 'x2')], type = "UK")
fine.grid$prob.delta.rando = pnorm(q = 0, mean = res$mean - 1, sd = res$sd)
##
# Threshold cutoff
fine.grid$cutof = 0
fine.grid$cutof[fine.grid$f1.norm <= 1 & fine.grid$f2.norm <= 1] = 1
mod.f1.paret = fill.sample.mod(GPar.data = data.paret, input.name = c('x1', 'x2'), output.name = 'f1.norm')
mod.f2.paret = fill.sample.mod(GPar.data = data.paret, input.name = c('x1', 'x2'), output.name = 'f2.norm')
mod.f1.adapt = fill.sample.mod(GPar.data = data.cutof, input.name = c('x1', 'x2'), output.name = 'f1.norm')
mod.f2.adapt = fill.sample.mod(GPar.data = data.cutof, input.name = c('x1', 'x2'), output.name = 'f2.norm')
mod.f1.rando = fill.sample.mod(GPar.data = data.rando, input.name = c('x1', 'x2'), output.name = 'f1.norm')
mod.f2.rando = fill.sample.mod(GPar.data = data.rando, input.name = c('x1', 'x2'), output.name = 'f2.norm')
res1 = predict(object = mod.f1.paret, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
res2 = predict(object = mod.f2.paret, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
fine.grid$prob.cutof.paret = pnorm(q = 0, mean = res1$mean - 1, sd = res1$sd) *
pnorm(q = 0, mean = res2$mean - 1, sd = res2$sd)
res1 = predict(object = mod.f1.adapt, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
res2 = predict(object = mod.f2.adapt, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
fine.grid$prob.cutof.adapt = pnorm(q = 0, mean = res1$mean - 1, sd = res1$sd) *
pnorm(q = 0, mean = res2$mean - 1, sd = res2$sd)
res1 = predict(object = mod.f1.rando, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
res2 = predict(object = mod.f2.rando, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
fine.grid$prob.cutof.rando = pnorm(q = 0, mean = res1$mean - 1, sd = res1$sd) *
pnorm(q = 0, mean = res2$mean - 1, sd = res2$sd)
##
# Radius-angle
fine.grid$radan = 0
fine.grid$radan[fine.grid$rad <= 1 & fine.grid$ang > 20] = 1
mod.rad.paret = fill.sample.mod(GPar.data = data.paret, input.name = c('x1', 'x2'), output.name = 'rad')
mod.ang.paret = fill.sample.mod(GPar.data = data.paret, input.name = c('x1', 'x2'), output.name = 'theta')
mod.rad.adapt = fill.sample.mod(GPar.data = data.radan, input.name = c('x1', 'x2'), output.name = 'rad')
mod.ang.adapt = fill.sample.mod(GPar.data = data.radan, input.name = c('x1', 'x2'), output.name = 'theta')
mod.rad.rando = fill.sample.mod(GPar.data = data.rando, input.name = c('x1', 'x2'), output.name = 'rad')
mod.ang.rando = fill.sample.mod(GPar.data = data.rando, input.name = c('x1', 'x2'), output.name = 'theta')
res1 = predict(object = mod.rad.paret, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
res2 = predict(object = mod.ang.paret, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
fine.grid$prob.radan.paret = pnorm(q = 0, mean = res1$mean - 1, sd = res1$sd) *
(1 - pnorm(q = 0, mean = res2$mean - 20, sd = res2$sd))
res1 = predict(object = mod.rad.adapt, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
res2 = predict(object = mod.ang.adapt, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
fine.grid$prob.radan.adapt = pnorm(q = 0, mean = res1$mean - 1, sd = res1$sd) *
(1 - pnorm(q = 0, mean = res2$mean - 20, sd = res2$sd))
res1 = predict(object = mod.rad.rando, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
res2 = predict(object = mod.ang.rando, newdata = data.frame(x1 = fine.grid$x1, x2 = fine.grid$x2), type = "UK")
fine.grid$prob.radan.rando = pnorm(q = 0, mean = res1$mean - 1, sd = res1$sd) *
(1 - pnorm(q = 0, mean = res2$mean - 20, sd = res2$sd))
sep = 0
ggplot() +
# Boundaries: +/- some separation from 0.5
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.delta.paret, color = 'delta'),
breaks = c(sep, -sep)+0.5) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.cutof.paret, color = 'cutof'),
breaks = c(sep, -sep)+0.5) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.radan.paret, color = 'radan'),
breaks = c(sep, -sep)+0.5) +
# Pareto frontier
geom_point(data = GPar.front, mapping = aes(x = x1, y = x2), color = 'black') +
geom_smooth(data = GPar.front, mapping = aes(x = x1, y = x2, color = 'Pareto'), level = 0.95, method = 'loess') +
labs(x = expression('x'[1]), y = expression('x'[2]), color = 'Acceptance Criteria', subtitle = 'Before Refinement') +
scale_color_manual(values = c('delta' = 'red', 'cutof' = 'skyblue2', 'radan' = 'green', 'Pareto' = 'black'),
labels = c('delta' = expression(delta*' < 1'), 'cutof' = expression('F'[1]^'*'*'< 1, F'[2]^'*'*'< 1'),
'radan' = expression('r < 1, '*theta*' > 18'^'o'),
'Pareto' = expression('Pareto Front')),
breaks = c('delta', 'cutof', 'radan', 'Pareto')) +
theme_classic() + theme(legend.position = c(0.85, 0.75)) +
scale_x_continuous(expand = c(0, 0), limits = c(0, 5)) + scale_y_continuous(expand = c(0, 0), limits = c(0, 5)) +
guides(colour = guide_legend(override.aes = list(fill = alpha('white', 1))))
ggplot() +
# Boundaries: +/- some separation from 0.5
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.delta.adapt, color = 'delta'),
breaks = c(sep, -sep)+0.5) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.cutof.adapt, color = 'cutof'),
breaks = c(sep, -sep)+0.5) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.radan.adapt, color = 'radan'),
breaks = c(sep, -sep)+0.5) +
# Pareto frontier
geom_point(data = GPar.front, mapping = aes(x = x1, y = x2), color = 'black') +
geom_smooth(data = GPar.front, mapping = aes(x = x1, y = x2, color = 'Pareto'), level = 0.95, method = 'loess') +
labs(x = expression('x'[1]), y = expression('x'[2]), color = 'Acceptance Criteria', subtitle = 'After Adaptive Sampling Refinement') +
scale_color_manual(values = c('delta' = 'red', 'cutof' = 'skyblue2', 'radan' = 'green', 'Pareto' = 'black'),
labels = c('delta' = expression(delta*' < 1'), 'cutof' = expression('F'[1]^'*'*'< 1, F'[2]^'*'*'< 1'),
'radan' = expression('r < 1, '*theta*' > 18'^'o'),
'Pareto' = expression('Pareto Front')),
breaks = c('delta', 'cutof', 'radan', 'Pareto')) +
theme_classic() + theme(legend.position = c(0.85, 0.75)) +
scale_x_continuous(expand = c(0, 0), limits = c(0, 5)) + scale_y_continuous(expand = c(0, 0), limits = c(0, 5)) +
guides(colour = guide_legend(override.aes = list(fill = alpha('white', 1))))
ggplot() +
# Boundaries: +/- some separation from 0.5
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.delta.rando, color = 'delta'),
breaks = c(sep, -sep)+0.5) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.cutof.rando, color = 'cutof'),
breaks = c(sep, -sep)+0.5) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.radan.rando, color = 'radan'),
breaks = c(sep, -sep)+0.5) +
# Pareto frontier
geom_point(data = GPar.front, mapping = aes(x = x1, y = x2), color = 'black') +
geom_smooth(data = GPar.front, mapping = aes(x = x1, y = x2, color = 'Pareto'), level = 0.95, method = 'loess') +
labs(x = expression('x'[1]), y = expression('x'[2]), color = 'Acceptance Criteria', subtitle = 'After Random Sampling') +
scale_color_manual(values = c('delta' = 'red', 'cutof' = 'skyblue2', 'radan' = 'green', 'Pareto' = 'black'),
labels = c('delta' = expression(delta*' < 1'), 'cutof' = expression('F'[1]^'*'*'< 1, F'[2]^'*'*'< 1'),
'radan' = expression('r < 1, '*theta*' > 18'^'o'),
'Pareto' = expression('Pareto Front')),
breaks = c('delta', 'cutof', 'radan', 'Pareto')) +
theme_classic() + theme(legend.position = c(0.85, 0.75)) +
scale_x_continuous(expand = c(0, 0), limits = c(0, 5)) + scale_y_continuous(expand = c(0, 0), limits = c(0, 5)) +
guides(colour = guide_legend(override.aes = list(fill = alpha('white', 1))))
# Plotting variables
sep = 0.05; trans = 0.25; ln.sz = 1.25
## Distance
ggplot() +
# Boundaries: +/- some separation from 0.5
geom_contour_filled(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.delta.paret,
fill = 'paret'),
breaks = c(sep, -sep)+0.5, linetype = 2, alpha = trans) +
geom_contour_filled(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.delta.adapt,
fill = 'adapt'),
breaks = c(sep, -sep)+0.5, linetype = 2, alpha = trans) +
geom_contour_filled(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.delta.rando,
fill = 'rando'),
breaks = c(sep, -sep)+0.5, linetype = 2, alpha = trans) +
# Actual boundaries
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.delta.paret, color = 'paret'),
breaks = 0.5, linetype = 1, size = ln.sz) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.delta.adapt, color = 'adapt'),
breaks = 0.5, linetype = 1, size = ln.sz) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.delta.rando, color = 'rando'),
breaks = 0.5, linetype = 1, size = ln.sz) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = delta, color = 'tru'),
breaks = 0.5, linetype = 1, size = ln.sz) +
scale_color_manual(breaks = c('tru', 'paret', 'adapt', 'rando'),
labels = c('tru' = 'True Boundary', 'paret' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'paret' = 'skyblue2', 'adapt' = 'red', 'rando' = 'green')) +
scale_fill_manual(breaks = c('tru', 'paret', 'adapt', 'rando'),
labels = c('tru' = 'True Boundary', 'paret' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'paret' = 'skyblue2', 'adapt' = 'red', 'rando' = 'green')) +
labs(x = expression('x'[1]), x = expression('x'[2]), subtitle = 'Criteria: Pareto Distance', color = '') +
theme_classic() + guides(fill = FALSE) + theme(legend.position = c(0.85, 0.85)) +
scale_x_continuous(limits = c(0, 5), expand = c(0, 0)) + scale_y_continuous(limits = c(0, 5), expand = c(0, 0))
## Threshold
ggplot() +
# Boundaries: +/- some separation from 0.5
geom_contour_filled(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.cutof.paret,
fill = 'paret'),
breaks = c(sep, -sep)+0.5, linetype = 2, alpha = trans) +
geom_contour_filled(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.cutof.adapt,
fill = 'adapt'),
breaks = c(sep, -sep)+0.5, linetype = 2, alpha = trans) +
geom_contour_filled(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.cutof.rando,
fill = 'rando'),
breaks = c(sep, -sep)+0.5, linetype = 2, alpha = trans) +
# Actual boundaries
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.cutof.paret, color = 'paret'),
breaks = 0.5, linetype = 1, size = ln.sz) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.cutof.adapt, color = 'adapt'),
breaks = 0.5, linetype = 1, size = ln.sz) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.cutof.rando, color = 'rando'),
breaks = 0.5, linetype = 1, size = ln.sz) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = cutof, color = 'tru'),
breaks = 0.5, linetype = 1, size = ln.sz) +
scale_color_manual(breaks = c('tru', 'paret', 'adapt', 'rando'),
labels = c('tru' = 'True Boundary', 'paret' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'paret' = 'skyblue2', 'adapt' = 'red', 'rando' = 'green')) +
scale_fill_manual(breaks = c('tru', 'paret', 'adapt', 'rando'),
labels = c('tru' = 'True Boundary', 'paret' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'paret' = 'skyblue2', 'adapt' = 'red', 'rando' = 'green')) +
labs(x = expression('x'[1]), x = expression('x'[2]), subtitle = 'Criteria: Objective function values', color = '') +
theme_classic() + guides(fill = FALSE) + theme(legend.position = c(0.85, 0.85)) +
scale_x_continuous(limits = c(0, 5), expand = c(0, 0)) + scale_y_continuous(limits = c(0, 5), expand = c(0, 0))
## Radius-angle
ggplot() +
# Boundaries: +/- some separation from 0.5
geom_contour_filled(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.radan.paret,
fill = 'paret'),
breaks = c(sep, -sep)+0.5, linetype = 2, alpha = trans) +
geom_contour_filled(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.radan.adapt,
fill = 'adapt'),
breaks = c(sep, -sep)+0.5, linetype = 2, alpha = trans) +
geom_contour_filled(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.radan.rando,
fill = 'rando'),
breaks = c(sep, -sep)+0.5, linetype = 2, alpha = trans) +
# Actual boundaries
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.radan.paret, color = 'paret'),
breaks = 0.5, linetype = 1, size = ln.sz) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.radan.adapt, color = 'adapt'),
breaks = 0.5, linetype = 1, size = ln.sz) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = prob.radan.rando, color = 'rando'),
breaks = 0.5, linetype = 1, size = ln.sz) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = radan, color = 'tru'),
breaks = 0.5, linetype = 1, size = ln.sz) +
scale_color_manual(breaks = c('tru', 'paret', 'adapt', 'rando'),
labels = c('tru' = 'True Boundary', 'paret' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'paret' = 'skyblue2', 'adapt' = 'red', 'rando' = 'green')) +
scale_fill_manual(breaks = c('tru', 'paret', 'adapt', 'rando'),
labels = c('tru' = 'True Boundary', 'paret' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'paret' = 'skyblue2', 'adapt' = 'red', 'rando' = 'green')) +
labs(x = expression('x'[1]), x = expression('x'[2]), color = '',
subtitle = expression('Criteria: Utopia point distance and f'[1]*' Priority > 80%')) +
theme_classic() + guides(fill = FALSE) + theme(legend.position = c(0.85, 0.85)) +
scale_x_continuous(limits = c(0, 5), expand = c(0, 0)) + scale_y_continuous(limits = c(0, 5), expand = c(0, 0))
NA
NA
The error rate in the GP model estimate should account for the built-in uncertainty in the model. The probability that the model gives the wrong result given the \((x_1, x_2)\) coordinates is the probability that it accepts the point when it should reject it, or vice versa. Since all \((x_1, x_2)\) are equally likely in this mathematical function, the error rate of the GP model is the average of these probabilities. There should be an evident decrease in the error rate between the pre- and post-refinement models.
The Monte Carlo uncertainty in the error rate is \(\sqrt{\frac{p(1-p)}{N}}\).
Since the fine grid calculation was already calculated, using that instead of a fully random assignment.
Also include the rate of type I (P[should reject | accepted]) and type II error (P[should accept | rejected]) error rates.
By Bayes rule, \(P[X|Y] = P[X, Y] P[X] / P[Y]\)
error.rate = data.frame(method = 'GP', source = rep(c('0paret', '1adapt', '2rando'), 3),
criteria = c(rep('Pareto Distance', 3),
rep('Objective Cutoff', 3),
rep('Utopia Distance + Priority', 3)))
error.rate$criteria = factor(error.rate$criteria, levels = c("Pareto Distance",
"Objective Cutoff", "Utopia Distance + Priority"))
# Total error rate
error.rate$rate =
c(mean(c(1 - filter(fine.grid, delta == 1)$prob.delta.paret, filter(fine.grid, delta == 0)$prob.delta.paret)),
mean(c(1 - filter(fine.grid, delta == 1)$prob.delta.adapt, filter(fine.grid, delta == 0)$prob.delta.adapt)),
mean(c(1 - filter(fine.grid, delta == 1)$prob.delta.rando, filter(fine.grid, delta == 0)$prob.delta.rando)),
mean(c(1 - filter(fine.grid, cutof == 1)$prob.cutof.paret, filter(fine.grid, cutof == 0)$prob.cutof.paret)),
mean(c(1 - filter(fine.grid, cutof == 1)$prob.cutof.adapt, filter(fine.grid, cutof == 0)$prob.cutof.adapt)),
mean(c(1 - filter(fine.grid, cutof == 1)$prob.cutof.rando, filter(fine.grid, cutof == 0)$prob.cutof.rando)),
mean(c(1 - filter(fine.grid, radan == 1)$prob.radan.paret, filter(fine.grid, radan == 0)$prob.radan.paret)),
mean(c(1 - filter(fine.grid, radan == 1)$prob.radan.adapt, filter(fine.grid, radan == 0)$prob.radan.adapt)),
mean(c(1 - filter(fine.grid, radan == 1)$prob.radan.rando, filter(fine.grid, radan == 0)$prob.radan.rando)))
error.rate$err = sqrt(error.rate$rate*(1-error.rate$rate)/nrow(fine.grid))
# Type 2 error = P[should accept | rejected]
error.rate$typ2 =
c(sum(1 - filter(fine.grid, delta == 1)$prob.delta.paret)*nrow(filter(fine.grid, delta == 1))/nrow(fine.grid)^2/
mean(1 - fine.grid$prob.delta.paret),
sum(1 - filter(fine.grid, delta == 1)$prob.delta.adapt)*nrow(filter(fine.grid, delta == 1))/nrow(fine.grid)^2/
mean(1 - fine.grid$prob.delta.adapt),
sum(1 - filter(fine.grid, delta == 1)$prob.delta.rando)*nrow(filter(fine.grid, delta == 1))/nrow(fine.grid)^2/
mean(1 - fine.grid$prob.delta.rando),
sum(1 - filter(fine.grid, cutof == 1)$prob.cutof.paret)*nrow(filter(fine.grid, cutof == 1))/nrow(fine.grid)^2/
mean(1 - fine.grid$prob.cutof.paret),
sum(1 - filter(fine.grid, cutof == 1)$prob.cutof.adapt)*nrow(filter(fine.grid, cutof == 1))/nrow(fine.grid)^2/
mean(1 - fine.grid$prob.cutof.adapt),
sum(1 - filter(fine.grid, cutof == 1)$prob.cutof.rando)*nrow(filter(fine.grid, cutof == 1))/nrow(fine.grid)^2/
mean(1 - fine.grid$prob.cutof.rando),
sum(1 - filter(fine.grid, radan == 1)$prob.radan.paret)*nrow(filter(fine.grid, radan == 1))/nrow(fine.grid)^2/
mean(1 - fine.grid$prob.radan.paret),
sum(1 - filter(fine.grid, radan == 1)$prob.radan.adapt)*nrow(filter(fine.grid, radan == 1))/nrow(fine.grid)^2/
mean(1 - fine.grid$prob.radan.adapt),
sum(1 - filter(fine.grid, radan == 1)$prob.radan.rando)*nrow(filter(fine.grid, radan == 1))/nrow(fine.grid)^2/
mean(1 - fine.grid$prob.radan.rando))
error.rate$typ2.err = sqrt(error.rate$typ2*(1-error.rate$typ2) / nrow(fine.grid))
# c(rep(nrow(filter(fine.grid, delta == 0)), 3),
# rep(nrow(filter(fine.grid, cutof == 0)), 3),
# rep(nrow(filter(fine.grid, radan == 0)), 3)))
# Type 1 error = P[should reject | accepted]
error.rate$typ1 =
c(sum(filter(fine.grid, delta == 0)$prob.delta.paret)*nrow(filter(fine.grid, delta == 0))/nrow(fine.grid)^2/
mean(fine.grid$prob.delta.paret),
sum(filter(fine.grid, delta == 0)$prob.delta.adapt)*nrow(filter(fine.grid, delta == 0))/nrow(fine.grid)^2/
mean(fine.grid$prob.delta.adapt),
sum(filter(fine.grid, delta == 0)$prob.delta.rando)*nrow(filter(fine.grid, delta == 0))/nrow(fine.grid)^2/
mean(fine.grid$prob.delta.rando),
sum(filter(fine.grid, cutof == 0)$prob.cutof.paret)*nrow(filter(fine.grid, cutof == 0))/nrow(fine.grid)^2/
mean(fine.grid$prob.cutof.paret),
sum(filter(fine.grid, cutof == 0)$prob.cutof.adapt)*nrow(filter(fine.grid, cutof == 0))/nrow(fine.grid)^2/
mean(fine.grid$prob.cutof.adapt),
sum(filter(fine.grid, cutof == 0)$prob.cutof.rando)*nrow(filter(fine.grid, cutof == 0))/nrow(fine.grid)^2/
mean(fine.grid$prob.cutof.rando),
sum(filter(fine.grid, radan == 0)$prob.radan.paret)*nrow(filter(fine.grid, radan == 0))/nrow(fine.grid)^2/
mean(fine.grid$prob.radan.paret),
sum(filter(fine.grid, radan == 0)$prob.radan.adapt)*nrow(filter(fine.grid, radan == 0))/nrow(fine.grid)^2/
mean(fine.grid$prob.radan.adapt),
sum(filter(fine.grid, radan == 0)$prob.radan.rando)*nrow(filter(fine.grid, radan == 0))/nrow(fine.grid)^2/
mean(fine.grid$prob.radan.rando))
error.rate$typ1.err = sqrt(error.rate$typ1*(1-error.rate$typ1) / nrow(fine.grid))
# c(rep(nrow(filter(fine.grid, delta == 0)), 3),
# rep(nrow(filter(fine.grid, cutof == 0)), 3),
# rep(nrow(filter(fine.grid, radan == 0)), 3)))
error.rate
g.tot = ggplot(error.rate) +
geom_col(mapping = aes(x = source, y = rate, fill = source)) +
geom_errorbar(mapping = aes(x = source, y = rate, ymin = rate - err, ymax = rate + err), width = 0.5) +
facet_wrap(~criteria, scales = 'free') +
scale_x_discrete(breaks = c(), labels = c()) +
labs(x = '', y = 'Total\nError Rate') +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
theme_classic() + scale_y_continuous(expand = expansion(mult = c(0, .1)))
g.ty1 = ggplot(error.rate) +
geom_col(mapping = aes(x = source, y = typ1, fill = source)) +
geom_errorbar(mapping = aes(x = source, y = typ1, ymin = typ1 - typ1.err, ymax = typ1 + typ1.err), width = 0.5) +
facet_wrap(~criteria, scales = 'free') +
scale_x_discrete(breaks = c(), labels = c()) +
labs(x = '', y = 'False Positive\nError Rate') +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
theme_classic() + scale_y_continuous(expand = expansion(mult = c(0, .1)))
g.ty2 = ggplot(error.rate) +
geom_col(mapping = aes(x = source, y = typ2, fill = source)) +
geom_errorbar(mapping = aes(x = source, y = typ2, ymin = typ2 - typ2.err, ymax = typ2 + typ2.err), width = 0.5) +
facet_wrap(~criteria, scales = 'free') +
scale_x_discrete(breaks = c(), labels = c()) +
labs(x = '', y = 'False Negative\nError Rate') +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
theme_classic() + scale_y_continuous(expand = expansion(mult = c(0, .1)))
(g.tot + guides(fill = FALSE)) / g.ty1 / (g.ty2 + guides(fill = FALSE))
write.csv(error.rate, 'ErrorRates.csv')
rm(g.tot, g.ty1, g.ty2)
A closer inspection of each GP-discovered boundary compared to the actual boundary based on fine-resolution function evaluation.
# Single variable conditionals compared to the expected conditionals by integration
Infer.plt = read.csv('../Ex_Quartic/Marginals_delta.csv')
ggplot() +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob - psd, color = cat), linetype = 2) +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob + psd, color = cat), linetype = 2) +
geom_path(data = filter(Infer.plt, cat != 'tru'), mapping = aes(x = x, y = prob, color = cat)) +
facet_wrap(var~., nrow = 2) +
scale_color_manual(labels = c('tru' = 'Expected Marginal', 'start' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'start' = 'skyblue2',
'adapt' = 'red', 'rando' = 'green'),
breaks = c('tru', 'start', 'adapt', 'rando')) +
guides(color = guide_legend(override.aes = list(linetype = c(2, 1, 1, 1)))) +
labs(x = '', y = 'Probability of Acceptance', subtitle = expression('Pareto Distance'), color = '') +
scale_y_continuous(expand = c(0, 0.05)) + scale_x_continuous(expand = c(0, 0)) +
theme_bw()
Infer.plt = read.csv('../Ex_Quartic/Marginals_cutof.csv')
ggplot() +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob - psd, color = cat), linetype = 2) +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob + psd, color = cat), linetype = 2) +
geom_path(data = filter(Infer.plt, cat != 'tru'), mapping = aes(x = x, y = prob, color = cat)) +
facet_wrap(var~., nrow = 2) +
scale_color_manual(labels = c('tru' = 'Expected Marginal', 'start' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'start' = 'skyblue2',
'adapt' = 'red', 'rando' = 'green'),
breaks = c('tru', 'start', 'adapt', 'rando')) +
guides(color = guide_legend(override.aes = list(linetype = c(2, 1, 1, 1)))) +
labs(x = '', y = 'Probability of Acceptance', subtitle = expression('Cutoff Thresholds'), color = '') +
scale_y_continuous(expand = c(0, 0.05)) + scale_x_continuous(expand = c(0, 0)) +
theme_bw()
Infer.plt = read.csv('../Ex_Quartic/Marginals_radan.csv')
ggplot() +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob - psd, color = cat), linetype = 2) +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob + psd, color = cat), linetype = 2) +
geom_path(data = filter(Infer.plt, cat != 'tru'), mapping = aes(x = x, y = prob, color = cat)) +
facet_wrap(var~., nrow = 2, labeller = label_parsed) +
scale_color_manual(labels = c('tru' = 'Expected Marginal', 'start' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'start' = 'skyblue2',
'adapt' = 'red', 'rando' = 'green'),
breaks = c('tru', 'start', 'adapt', 'rando')) +
guides(color = guide_legend(override.aes = list(linetype = c(2, 1, 1, 1)))) +
labs(x = '', y = 'Probability of Acceptance', subtitle = expression('Utopia Distance + Priority'), color = '') +
scale_y_continuous(expand = c(0, 0.05)) + scale_x_continuous(expand = c(0, 0)) +
theme_bw()
NA
NA
The method should use only the data from the initial Pareto front calculation. Additional samples may be necessary to refine the results, but there is no established method for this refinement as of yet.
SVM methods require existing labels for the points; in this case, the labels are whether or not the point meets the same three selection criteria as tested with the new GP method.
Sigmoidal kernels always give very bad fits, so ignoring those possibilities in the error calculation.
# Create a function to output the plot of all 4 models compared to the real result
SVM.input = function(fine.input, mod.rad, mod.lin, mod.pol, mod.sig, tit){
# Fine input has the variables x1, x2, and cat
res = predict(object = mod.lin, newdata = fine.input[, c('x1', 'x2')])
fine.input$res = res
g.lin = ggplot(fine.input) +
geom_point(data = fine.input, mapping = aes(x = x1, y = x2, color = res)) +
geom_contour(data = fine.input, mapping = aes(x = x1, y = x2, z = cat), breaks = c(0.5)) +
labs(subtitle = 'Linear', x = expression('x'[1]), y = expression('x'[2])) +
guides(color = FALSE)
res = predict(object = mod.rad, newdata = fine.input[,c('x1', 'x2')])
fine.input$res = res
g.rad = ggplot() +
geom_point(data = fine.input, mapping = aes(x = x1, y = x2, color = res)) +
geom_contour(data = fine.input, mapping = aes(x = x1, y = x2, z = cat), breaks = c(0.5)) +
labs(subtitle = 'Radial', x = expression('x'[1]), y = expression('x'[2]), title = tit) + guides(color = FALSE)
res = predict(object = mod.pol, newdata = fine.input[,c('x1', 'x2')])
fine.input$res = res
g.pol = ggplot() +
geom_point(data = fine.input, mapping = aes(x = x1, y = x2, color = res)) +
geom_contour(data = fine.input, mapping = aes(x = x1, y = x2, z = cat), breaks = c(0.5)) +
labs(subtitle = 'Polynomial', x = expression('x'[1]), y = expression('x'[2])) + guides(color = FALSE)
res = predict(object = mod.sig, newdata = fine.input[,c('x1', 'x2')])
fine.input$res = res
g.sig = ggplot() +
geom_point(data = fine.input, mapping = aes(x = x1, y = x2, color = res)) +
geom_contour(data = fine.input, mapping = aes(x = x1, y = x2, z = cat), breaks = c(0.5)) +
labs(subtitle = 'Sigmoid', x = expression('x'[1]), y = expression('x'[2])) + guides(color = FALSE)
return((g.rad + g.lin) / (g.pol + g.sig))
}
SVM.mod.all = function(input.data){
# Given input data, find all 4 default models and return as a list
fit.rad = svm(as.factor(cat) ~ x1*x2, data = input.data[,c('x1', 'x2', 'cat')],
scale = FALSE, kernel = "radial", cost = 5)
fit.lin = svm(as.factor(cat) ~ x1*x2, data = input.data[,c('x1', 'x2', 'cat')],
scale = FALSE, kernel = "linear", cost = 5)
# Polynomial: increase the kernel degree due to complexity of the boundary; maximum tested that still achieved convergence
fit.pol = svm(as.factor(cat) ~ x1*x2, data = input.data[,c('x1', 'x2', 'cat')],
scale = FALSE, kernel = "polynomial", cost = 5, degree = 3.5)
# Increease the coefficient
fit.sig = svm(as.factor(cat) ~ x1*x2, data = input.data[,c('x1', 'x2', 'cat')],
scale = FALSE, kernel = "sigmoid", cost = 5)
return(list(rad = fit.rad, lin = fit.lin, pol = fit.pol, sig = fit.sig))
}
SVM.err = function(fine.input, mod.list){
# Fine grid has the real result; mod.list is the 4 different SVM kernels
fine.input$rad = predict(object = mod.list$rad, newdata = fine.input[,c('x1', 'x2')])
fine.input$lin = predict(object = mod.list$lin, newdata = fine.input[,c('x1', 'x2')])
fine.input$pol = predict(object = mod.list$pol, newdata = fine.input[,c('x1', 'x2')])
# fine.input$sig = predict(object = mod.list$sig, newdata = fine.input[,c('x1', 'x2')])
# Error rate
rad = nrow(filter(fine.input, cat == 0, rad == 1)) + nrow(filter(fine.input, cat == 1, rad == 0))
lin = nrow(filter(fine.input, cat == 0, lin == 1)) + nrow(filter(fine.input, cat == 1, lin == 0))
pol = nrow(filter(fine.input, cat == 0, pol == 1)) + nrow(filter(fine.input, cat == 1, pol == 0))
# sig = nrow(filter(fine.input, cat == 0, sig == 1)) + nrow(filter(fine.input, cat == 1, sig == 0))
# rate = c(rad, lin, pol, sig)/nrow(fine.input)
rate = c(rad, lin, pol)/nrow(fine.input)
# Uncertainty
err = sqrt(rate*(1-rate)/nrow(fine.input))
# Type I error = P[should reject | accepted] = P[both] P[should reject] / P[accepted]
typ1 =
c(nrow(filter(fine.input, cat == 0, rad == 1))*nrow(filter(fine.input, cat == 0))/nrow(fine.input)/
nrow(filter(fine.input, rad == 1)),
nrow(filter(fine.input, cat == 0, lin == 1))*nrow(filter(fine.input, cat == 0))/nrow(fine.input)/
nrow(filter(fine.input, lin == 1)),
nrow(filter(fine.input, cat == 0, pol == 1))*nrow(filter(fine.input, cat == 0))/nrow(fine.input)/
nrow(filter(fine.input, pol == 1)))#,
# nrow(filter(fine.input, cat == 0, sig == 1))*nrow(filter(fine.input, cat == 0))/nrow(fine.input)/
# nrow(filter(fine.input, sig == 1)) )
typ1.err = sqrt(typ1*(1-typ1) / nrow(fine.input))
# Type II error = P[should accept | rejected] = P[both] P[should accept] / P[rejected]
typ2 =
c(nrow(filter(fine.input, cat == 1, rad == 0))*nrow(filter(fine.input, cat == 1))/nrow(fine.input)/
nrow(filter(fine.input, rad == 0)),
nrow(filter(fine.input, cat == 1, lin == 0))*nrow(filter(fine.input, cat == 1))/nrow(fine.input)/
nrow(filter(fine.input, lin == 0)),
nrow(filter(fine.input, cat == 1, pol == 0))*nrow(filter(fine.input, cat == 1))/nrow(fine.input)/
nrow(filter(fine.input, pol == 0)))#,
# nrow(filter(fine.input, cat == 1, sig == 0))*nrow(filter(fine.input, cat == 1))/nrow(fine.input)/
# nrow(filter(fine.input, sig == 0)) )
typ2.err = sqrt(typ2*(1-typ2) / nrow(fine.input))
# return(data.frame(method = c('SVM-rad', 'SVM-lin', 'SVM-pol', 'SVM-sig'), rate, err, typ1, typ1.err, typ2, typ2.err))
return(data.frame(method = c('SVM-rad', 'SVM-lin', 'SVM-pol'), rate, err, typ1, typ1.err, typ2, typ2.err))
}
# For the radial and polynomial kernels, calculate the marginals and Shapley values for comparison
marginal = function(mod){
x.rng = seq(from = 0, to = 5, length.out = 50)
svm.margin = data.frame()
nsamp = 2000
for(x in x.rng){
# x1 marginal
temp.frame = data.frame(x1 = x, x2 = runif(n = nsamp, min = 0, max = 5))
res = as.numeric(predict(object = mod$pol, newdata = temp.frame))-1
svm.margin = rbind(svm.margin, data.frame(x = x, var = 'x1', prob = sum(res)/length(res), method = 'SVM-pol'))
res = as.numeric(predict(object = mod$rad, newdata = temp.frame))-1
svm.margin = rbind(svm.margin, data.frame(x = x, var = 'x1', prob = sum(res)/length(res), method = 'SVM-rad'))
# x1 marginal
temp.frame = data.frame(x2 = x, x1 = runif(n = nsamp, min = 0, max = 5))
res = as.numeric(predict(object = mod$pol, newdata = temp.frame))-1
svm.margin = rbind(svm.margin, data.frame(x = x, var = 'x2', prob = sum(res)/length(res), method = 'SVM-pol'))
res = as.numeric(predict(object = mod$rad, newdata = temp.frame))-1
svm.margin = rbind(svm.margin, data.frame(x = x, var = 'x2', prob = sum(res)/length(res), method = 'SVM-rad'))
}
svm.margin$psd = sqrt(svm.margin$prob*(1 - svm.margin$prob)/nsamp)
return(svm.margin)
}
SVM.shap = function(mod){
# Strumbelj et al. (2014) Monte Carlo estimate.
# Since this is a classification setting, only a difference of 0 or 1 is possible
# i.e. a difference of -1 is the same as a difference of 1, as it is simply a misclassification
nsamp = 1500*50 # Same number of samples as other methods for consistency
x0 = data.frame(x1 = runif(n = nsamp, min = 0, max = 5),
x2 = runif(n = nsamp, min = 0, max = 5))
z1 = data.frame(x1 = x0$x1,
x2 = runif(n = nsamp, min = 0, max = 5))
z2 = data.frame(x1 = runif(n = nsamp, min = 0, max = 5),
x2 = x0$x2)
# Apply functions to all
res.pol = as.numeric(predict(object = mod$pol, newdata = x0))
res.rad = as.numeric(predict(object = mod$rad, newdata = x0))
x0$pol = res.pol; x0$rad = res.rad
res.pol = as.numeric(predict(object = mod$pol, newdata = z1))
res.rad = as.numeric(predict(object = mod$rad, newdata = z1))
z1$pol = res.pol; z1$rad = res.rad
res.pol = as.numeric(predict(object = mod$pol, newdata = z2))
res.rad = as.numeric(predict(object = mod$rad, newdata = z2))
z2$pol = res.pol; z2$rad = res.rad
# Calculate mean and standard error of the differences for Shapley values
shap = c(mean((x0$pol - z1$pol)), mean((x0$pol - z2$pol)),
mean((x0$rad - z1$rad)), mean((x0$rad - z2$rad)))
r.shap = c(shap[1:2]/max(shap[1:2]), shap[3:4]/max(shap[3:4]))
# For standard error, calculate standard error of the mean with n = 1500
# to be consistent with the other standard errors
pol.x1m = apply(matrix(x0$pol - z1$pol, nrow = 1500), 2, mean)
rad.x1m = apply(matrix(x0$rad - z1$rad, nrow = 1500), 2, mean)
pol.x2m = apply(matrix(x0$pol - z2$pol, nrow = 1500), 2, mean)
rad.x2m = apply(matrix(x0$rad - z2$rad, nrow = 1500), 2, mean)
shap.err = c(sd(pol.x1m), sd(pol.x2m),
sd(rad.x1m), sd(rad.x2m))
return(data.frame(import = shap, var = c('x1', 'x2'),
sd = shap.err, r.import = r.shap,
method = c('Shapley-pol', 'Shapley-pol',
'Shapley-rad', 'Shapley-rad')))
}
data.paret = read.csv('../Ex_Quartic/GPar_all_start.csv')
data.delta = read.csv('../Ex_Quartic/GPar_Accept_Delta1.csv')
data.rando = read.csv(file = 'GPar_Random.csv')
# Define acceptance
data.paret$cat = 0; data.delta$cat = 0; data.rando$cat = 0; fine.grid$cat = 0
data.paret$cat[data.paret$dist <= 1] = 1
data.delta$cat[data.delta$dist <= 1] = 1
data.rando$cat[data.rando$dist <= 1] = 1
fine.grid$cat[fine.grid$dist <= 1] = 1
# Plots
mod.paret = SVM.mod.all(input.data = data.paret)
SVM.input(fine.input = fine.grid[,c('x1', 'x2', 'cat')],
mod.rad = mod.paret$rad, mod.lin = mod.paret$lin,
mod.pol = mod.paret$pol, mod.sig = mod.paret$sig,
tit = 'Pareto Distance, Starting Dataset')
mod.adapt = SVM.mod.all(input.data = data.delta)
# Store this model for importance ranking later
mod.adapt.delta = mod.adapt
SVM.input(fine.input = fine.grid[,c('x1', 'x2', 'cat')],
mod.rad = mod.adapt$rad, mod.lin = mod.adapt$lin,
mod.pol = mod.adapt$pol, mod.sig = mod.adapt$sig,
tit = 'Pareto Distance, + Adaptive Sampling')
mod.rando = SVM.mod.all(input.data = data.rando)
SVM.input(fine.input = fine.grid[,c('x1', 'x2', 'cat')],
mod.rad = mod.rando$rad, mod.lin = mod.rando$lin,
mod.pol = mod.rando$pol, mod.sig = mod.rando$sig,
tit = 'Pareto Distance, + Random Sampling')
err.paret = SVM.err(fine.input = fine.grid, mod.list = mod.paret)
err.paret$source = '0paret'; err.paret$criteria = 'Pareto Distance'
err.adapt = SVM.err(fine.input = fine.grid, mod.list = mod.adapt)
err.adapt$source = '1adapt'; err.adapt$criteria = 'Pareto Distance'
err.rando = SVM.err(fine.input = fine.grid, mod.list = mod.rando)
err.rando$source = '2rando'; err.rando$criteria = 'Pareto Distance'
err.svm.dist = rbind(err.paret, err.adapt, err.rando)
g.rate = ggplot(err.svm.dist) +
geom_col(mapping = aes(x = source, y = rate, fill = source)) +
geom_errorbar(mapping = aes(x = source, ymin = rate - err, ymax = rate + err), width = 0.5) +
facet_wrap(~method, nrow = 1) +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + scale_x_discrete(breaks = '') +
theme_bw() + labs(x = '', y = 'Total\nError Rate', subtitle = 'Pareto Distance')
g.typ1 = ggplot(err.svm.dist) +
geom_col(mapping = aes(x = source, y = typ1, fill = source)) +
geom_errorbar(mapping = aes(x = source, ymin = typ1 - typ1.err, ymax = typ1 + typ1.err), width = 0.5) +
facet_wrap(~method, nrow = 1) +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + scale_x_discrete(breaks = '') +
theme_bw() + labs(x = '', y = 'False Positive\nError Rate')
g.typ2 = ggplot(err.svm.dist) +
geom_col(mapping = aes(x = source, y = typ2, fill = source)) +
geom_errorbar(mapping = aes(x = source, ymin = typ2 - typ2.err, ymax = typ2 + typ2.err), width = 0.5) +
facet_wrap(~method, nrow = 1) +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + scale_x_discrete(breaks = '') +
theme_bw() + labs(x = '', y = 'False Negative\nError Rate')
(g.rate + guides(fill = FALSE)) / g.typ1 / (g.typ2 + guides(fill = FALSE))
write.csv(err.svm.dist, 'SVM_error_delta.csv')
rm(err.paret, err.adapt, err.rando)
Graphically, the selection of the kernel form has a large impact on the quality of the result. The polynomial and radial kernels appear to roughly match the shape of the normalized distance = 1 criterion; the linear and sigmoidal kernels do not appear to match anything. Further tuning of the polynomial degree or coefficient could improve the results, but without a training-testing set delineation, there is no way to validate the results.
svm.margin.paret = marginal(mod = mod.paret); svm.margin.paret$cat = 'start'
svm.margin.adapt = marginal(mod = mod.adapt); svm.margin.adapt$cat = 'adapt'
svm.margin.rando = marginal(mod = mod.rando); svm.margin.rando$cat = 'rando'
svm.margin = rbind(svm.margin.paret, svm.margin.adapt, svm.margin.rando)
rm(svm.margin.paret, svm.margin.adapt, svm.margin.rando)
write.csv(svm.margin, 'Margin_SVM_dist.csv')
Infer.plt = read.csv('../Ex_Quartic/Marginals_delta.csv')
Infer.plt = Infer.plt[,!names(Infer.plt) %in% c('X')]
# names(Infer.plt)
Infer.plt$method = 'GP'
ggplot() +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob - psd, color = cat), linetype = 3) +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob + psd, color = cat), linetype = 3) +
geom_path(data = svm.margin, mapping = aes(x = x, y = prob, linetype = method, color = cat)) +
facet_wrap(~var, nrow = 2) + theme_bw() +
labs(x = 'Input Value', y = 'Conditional Probability', linetype = 'SVM Kernel', color = '') +
scale_color_manual(labels = c('tru' = 'Expected Marginal', 'start' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'start' = 'skyblue2',
'adapt' = 'red', 'rando' = 'green'),
breaks = c('tru', 'start', 'adapt', 'rando')) +
guides(color = guide_legend(override.aes = list(linetype = c(3, 1, 1, 1)))) +
scale_linetype_discrete(labels = c('SVM-pol' = 'Polynomial', 'rad' = 'SVM-Radial')) +
scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(limits = c(0, 1))
ggplot() +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob - psd, color = cat), linetype = 3) +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob + psd, color = cat), linetype = 3) +
geom_path(data = filter(Infer.plt, cat != 'tru'), mapping = aes(x = x, y = prob, color = cat)) +
facet_wrap(var~., nrow = 2) +
scale_color_manual(labels = c('tru' = 'Expected Marginal', 'start' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'start' = 'skyblue2',
'adapt' = 'red', 'rando' = 'green'),
breaks = c('tru', 'start', 'adapt', 'rando')) +
guides(color = guide_legend(override.aes = list(linetype = c(3, 1, 1, 1)))) +
labs(x = '', y = 'Probability of Acceptance', subtitle = expression('Pareto Distance'), color = '') +
scale_y_continuous(expand = c(0, 0.05)) + scale_x_continuous(expand = c(0, 0)) +
theme_bw()
# Calculate coefficients of determination: easier comparison
coefdet = data.frame(method = c(rep('SVM-pol', 3), rep('SVM-rad', 3), rep('GP', 3)),
dataset = c('0start', '1adapt', '2rando'),
coef = c(cor(x = filter(svm.margin, method == 'SVM-pol', cat == 'start')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-pol', cat == 'adapt')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-pol', cat == 'rando')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-rad', cat == 'start')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-rad', cat == 'adapt')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-rad', cat == 'rando')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(Infer.plt, cat == 'start')$prob,
y = c(filter(Infer.plt, cat == 'tru', var == 'x1')$prob, filter(Infer.plt, cat == 'tru', var == 'x2')$prob),
method = 'pearson')^2,
cor(x = filter(Infer.plt, cat == 'adapt')$prob,
y = c(filter(Infer.plt, cat == 'tru', var == 'x1')$prob, filter(Infer.plt, cat == 'tru', var == 'x2')$prob),
method = 'pearson')^2,
cor(x = filter(Infer.plt, cat == 'rando')$prob,
y = c(filter(Infer.plt, cat == 'tru', var == 'x1')$prob, filter(Infer.plt, cat == 'tru', var == 'x2')$prob),
method = 'pearson')^2))
coefdet
write.csv(coefdet, 'Marginals_CoefDet_delta.csv')
ggplot(coefdet) +
geom_col(mapping = aes(x = dataset, fill = dataset, y = coef))+
facet_grid(~method) +
labs(x = '', y = '1-Variable Marginal Coefficient of Determination', subtitle = 'Pareto Distance') +
scale_x_discrete(breaks = c()) +
scale_fill_manual(labels = c('0start' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0start' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05)))
Single variable marginals are ok visually, but calculating the coefficient of determination shows that the GP method is more consistent regardless of dataset and it demonstrates the expected improvement with increased sampling near the boundary. The SVM methods overall cannot capture the boundary very well, leading to classification errors.
data.paret = read.csv('../Ex_Quartic/GPar_all_start.csv')
data.cutof = read.csv('../Ex_Quartic/GPar_Accept_Threshold.csv')
data.rando = read.csv(file = 'GPar_Random.csv')
# Define acceptance
data.paret$cat = 0; data.cutof$cat = 0; data.rando$cat = 0; fine.grid$cat = 0
data.paret$cat[data.paret$f1.norm <= 1 & data.paret$f2.norm <= 1] = 1
data.cutof$cat[data.cutof$f1.norm <= 1 & data.cutof$f2.norm <= 1] = 1
data.rando$cat[data.rando$f1.norm <= 1 & data.rando$f2.norm <= 1] = 1
fine.grid$cat[fine.grid$f1.norm <= 1 & fine.grid$f2.norm <= 1] = 1
# Plots
mod.paret = SVM.mod.all(input.data = data.paret)
SVM.input(fine.input = fine.grid[,c('x1', 'x2', 'cat')],
mod.rad = mod.paret$rad, mod.lin = mod.paret$lin,
mod.pol = mod.paret$pol, mod.sig = mod.paret$sig,
tit = 'Threshold Cutoff, Starting Dataset')
mod.adapt = SVM.mod.all(input.data = data.cutof)
mod.adapt.cutof = mod.adapt
SVM.input(fine.input = fine.grid[,c('x1', 'x2', 'cat')],
mod.rad = mod.adapt$rad, mod.lin = mod.adapt$lin,
mod.pol = mod.adapt$pol, mod.sig = mod.adapt$sig,
tit = 'Threshold Cutoff, + Adaptive Sampling')
mod.rando = SVM.mod.all(input.data = data.rando)
SVM.input(fine.input = fine.grid[,c('x1', 'x2', 'cat')],
mod.rad = mod.rando$rad, mod.lin = mod.rando$lin,
mod.pol = mod.rando$pol, mod.sig = mod.rando$sig,
tit = 'Threshold Cutoff, + Random Sampling')
# Error rates
err.paret = SVM.err(fine.input = fine.grid, mod.list = mod.paret)
err.paret$source = '0paret'; err.paret$criteria = 'Treshold Cutoff'
err.adapt = SVM.err(fine.input = fine.grid, mod.list = mod.adapt)
err.adapt$source = '1adapt'; err.adapt$criteria = 'Treshold Cutoff'
err.rando = SVM.err(fine.input = fine.grid, mod.list = mod.rando)
err.rando$source = '2rando'; err.rando$criteria = 'Treshold Cutoff'
err.svm.cutof = rbind(err.paret, err.adapt, err.rando)
g.rate = ggplot(err.svm.cutof) +
geom_col(mapping = aes(x = source, y = rate, fill = source)) +
geom_errorbar(mapping = aes(x = source, ymin = rate - err, ymax = rate + err), width = 0.5) +
facet_wrap(~method, nrow = 1) +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + scale_x_discrete(breaks = '') +
theme_bw() + labs(x = '', y = 'Total\nError Rate', subtitle = 'Cutoff Threshold')
g.typ1 = ggplot(err.svm.cutof) +
geom_col(mapping = aes(x = source, y = typ1, fill = source)) +
geom_errorbar(mapping = aes(x = source, ymin = typ1 - typ1.err, ymax = typ1 + typ1.err), width = 0.5) +
facet_wrap(~method, nrow = 1) +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + scale_x_discrete(breaks = '') +
theme_bw() + labs(x = '', y = 'False Positive\nError Rate')
g.typ2 = ggplot(err.svm.cutof) +
geom_col(mapping = aes(x = source, y = typ2, fill = source)) +
geom_errorbar(mapping = aes(x = source, ymin = typ2 - typ2.err, ymax = typ2 + typ2.err), width = 0.5) +
facet_wrap(~method, nrow = 1) +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + scale_x_discrete(breaks = '') +
theme_bw() + labs(x = '', y = 'False Negative\nError Rate')
(g.rate + guides(fill = FALSE)) / g.typ1 / (g.typ2 + guides(fill = FALSE))
write.csv(err.svm.cutof, 'SVM_error_cutof.csv')
rm(err.paret, err.adapt, err.rando, g.rate, g.typ1, g.typ2)
svm.margin.paret = marginal(mod = mod.paret); svm.margin.paret$cat = 'start'
svm.margin.adapt = marginal(mod = mod.adapt); svm.margin.adapt$cat = 'adapt'
svm.margin.rando = marginal(mod = mod.rando); svm.margin.rando$cat = 'rando'
svm.margin = rbind(svm.margin.paret, svm.margin.adapt, svm.margin.rando)
rm(svm.margin.paret, svm.margin.adapt, svm.margin.rando)
write.csv(svm.margin, 'Margin_SVM_cutof.csv')
Infer.plt = read.csv('../Ex_Quartic/Marginals_cutof.csv')
Infer.plt = Infer.plt[,!names(Infer.plt) %in% c('X')]
# names(Infer.plt)
Infer.plt$method = 'GP'
ggplot() +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob - psd, color = cat), linetype = 3) +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob + psd, color = cat), linetype = 3) +
geom_path(data = svm.margin, mapping = aes(x = x, y = prob, linetype = method, color = cat)) +
facet_wrap(~var, nrow = 2) + theme_bw() +
labs(x = 'Input Value', y = 'Conditional Probability', linetype = 'SVM Kernel', color = '', subtitle = 'Cutoff Threshold') +
scale_color_manual(labels = c('tru' = 'Expected Marginal', 'start' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'start' = 'skyblue2',
'adapt' = 'red', 'rando' = 'green'),
breaks = c('tru', 'start', 'adapt', 'rando')) +
guides(color = guide_legend(override.aes = list(linetype = c(3, 1, 1, 1)))) +
scale_linetype_discrete(labels = c('SVM-pol' = 'Polynomial', 'rad' = 'SVM-Radial')) +
scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(limits = c(0, 1))
ggplot() +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob - psd, color = cat), linetype = 3) +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob + psd, color = cat), linetype = 3) +
geom_path(data = filter(Infer.plt, cat != 'tru'), mapping = aes(x = x, y = prob, color = cat)) +
facet_wrap(var~., nrow = 2) +
scale_color_manual(labels = c('tru' = 'Expected Marginal', 'start' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'start' = 'skyblue2',
'adapt' = 'red', 'rando' = 'green'),
breaks = c('tru', 'start', 'adapt', 'rando')) +
guides(color = guide_legend(override.aes = list(linetype = c(3, 1, 1, 1)))) +
labs(x = '', y = 'Probability of Acceptance', subtitle = 'Cutoff Threshold', color = '') +
scale_y_continuous(expand = c(0, 0.05)) + scale_x_continuous(expand = c(0, 0)) +
theme_bw()
# Calculate coefficients of determination: easier comparison
coefdet = data.frame(method = c(rep('SVM-pol', 3), rep('SVM-rad', 3), rep('GP', 3)),
dataset = c('0start', '1adapt', '2rando'),
coef = c(cor(x = filter(svm.margin, method == 'SVM-pol', cat == 'start')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-pol', cat == 'adapt')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-pol', cat == 'rando')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-rad', cat == 'start')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-rad', cat == 'adapt')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-rad', cat == 'rando')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(Infer.plt, cat == 'start')$prob,
y = c(filter(Infer.plt, cat == 'tru', var == 'x1')$prob, filter(Infer.plt, cat == 'tru', var == 'x2')$prob),
method = 'pearson')^2,
cor(x = filter(Infer.plt, cat == 'adapt')$prob,
y = c(filter(Infer.plt, cat == 'tru', var == 'x1')$prob, filter(Infer.plt, cat == 'tru', var == 'x2')$prob),
method = 'pearson')^2,
cor(x = filter(Infer.plt, cat == 'rando')$prob,
y = c(filter(Infer.plt, cat == 'tru', var == 'x1')$prob, filter(Infer.plt, cat == 'tru', var == 'x2')$prob),
method = 'pearson')^2))
coefdet
write.csv(coefdet, 'Marginals_CoefDet_cutof.csv')
ggplot(coefdet) +
geom_col(mapping = aes(x = dataset, fill = dataset, y = coef))+
facet_grid(~method) +
labs(x = '', y = '1-Variable Marginal Coefficient of Determination', subtitle = 'Cutoff Threshold') +
scale_x_discrete(breaks = c()) +
scale_fill_manual(labels = c('0start' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0start' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05)))
Repeat for radius-angle cutoff
data.paret = read.csv(file = '../Ex_Quartic/GPar_all_start.csv')
data.radan = read.csv('../Ex_Quartic/GPar_Accept_Radius.csv')
data.rando = read.csv(file = 'GPar_Random.csv')
# Define acceptance
data.paret$cat = 0; data.radan$cat = 0; data.rando$cat = 0; fine.grid$cat = 0
data.paret$cat[sqrt(data.paret$f1.norm^2 + data.paret$f2.norm^2) <= 1 & data.paret$theta > 20] = 1
data.radan$cat[sqrt(data.radan$f1.norm^2 + data.radan$f2.norm^2) <= 1 & data.radan$theta > 20] = 1
data.rando$cat[sqrt(data.rando$f1.norm^2 + data.rando$f2.norm^2) <= 1 & data.rando$theta > 20] = 1
fine.grid$cat[sqrt(fine.grid$f1.norm^2 + fine.grid$f2.norm^2) <= 1 & fine.grid$ang > 20] = 1
# Plots
mod.paret = SVM.mod.all(input.data = data.paret)
SVM.input(fine.input = fine.grid[,c('x1', 'x2', 'cat')],
mod.rad = mod.paret$rad, mod.lin = mod.paret$lin,
mod.pol = mod.paret$pol, mod.sig = mod.paret$sig,
tit = 'Utopia Distance + Priority, Starting Dataset')
mod.adapt = SVM.mod.all(input.data = data.radan)
mod.adapt.radan = mod.adapt
SVM.input(fine.input = fine.grid[,c('x1', 'x2', 'cat')],
mod.rad = mod.adapt$rad, mod.lin = mod.adapt$lin,
mod.pol = mod.adapt$pol, mod.sig = mod.adapt$sig,
tit = 'Utopia Distance + Priority, + Adaptive Sampling')
mod.rando = SVM.mod.all(input.data = data.rando)
SVM.input(fine.input = fine.grid[,c('x1', 'x2', 'cat')],
mod.rad = mod.rando$rad, mod.lin = mod.rando$lin,
mod.pol = mod.rando$pol, mod.sig = mod.rando$sig,
tit = 'Utopia Distance + Priority, + Random Sampling')
# Error rates
err.paret = SVM.err(fine.input = fine.grid, mod.list = mod.paret)
err.paret$source = '0paret'; err.paret$criteria = 'Treshold Cutoff'
err.adapt = SVM.err(fine.input = fine.grid, mod.list = mod.adapt)
err.adapt$source = '1adapt'; err.adapt$criteria = 'Treshold Cutoff'
err.rando = SVM.err(fine.input = fine.grid, mod.list = mod.rando)
err.rando$source = '2rando'; err.rando$criteria = 'Treshold Cutoff'
err.svm.radan = rbind(err.paret, err.adapt, err.rando)
g.rate = ggplot(err.svm.radan) +
geom_col(mapping = aes(x = source, y = rate, fill = source)) +
geom_errorbar(mapping = aes(x = source, ymin = rate - err, ymax = rate + err), width = 0.5) +
facet_wrap(~method, nrow = 1) +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + scale_x_discrete(breaks = '') +
theme_bw() + labs(x = '', y = 'Total\nError Rate', subtitle = 'Utopia Distance + Priority')
g.typ1 = ggplot(err.svm.radan) +
geom_col(mapping = aes(x = source, y = typ1, fill = source)) +
geom_errorbar(mapping = aes(x = source, ymin = typ1 - typ1.err, ymax = typ1 + typ1.err), width = 0.5) +
facet_wrap(~method, nrow = 1) +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + scale_x_discrete(breaks = '') +
theme_bw() + labs(x = '', y = 'False Positive\nError Rate')
g.typ2 = ggplot(err.svm.radan) +
geom_col(mapping = aes(x = source, y = typ2, fill = source)) +
geom_errorbar(mapping = aes(x = source, ymin = typ2 - typ2.err, ymax = typ2 + typ2.err), width = 0.5) +
facet_wrap(~method, nrow = 1) +
scale_fill_manual(labels = c('0paret' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0paret' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + scale_x_discrete(breaks = '') +
theme_bw() + labs(x = '', y = 'False Negative\nError Rate')
(g.rate + guides(fill = FALSE)) / g.typ1 / (g.typ2 + guides(fill = FALSE))
write.csv(err.svm.cutof, 'SVM_error_radan.csv')
rm(err.paret, err.adapt, err.rando, g.rate, g.typ1, g.typ2)
svm.margin.paret = marginal(mod = mod.paret); svm.margin.paret$cat = 'start'
svm.margin.adapt = marginal(mod = mod.adapt); svm.margin.adapt$cat = 'adapt'
svm.margin.rando = marginal(mod = mod.rando); svm.margin.rando$cat = 'rando'
svm.margin = rbind(svm.margin.paret, svm.margin.adapt, svm.margin.rando)
rm(svm.margin.paret, svm.margin.adapt, svm.margin.rando)
write.csv(svm.margin, 'Margin_SVM_radan.csv')
Infer.plt = read.csv('../Ex_Quartic/Marginals_radan.csv')
Infer.plt = Infer.plt[,!names(Infer.plt) %in% c('X')]
# names(Infer.plt)
Infer.plt$method = 'GP'
ggplot() +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob - psd, color = cat), linetype = 3) +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob + psd, color = cat), linetype = 3) +
geom_path(data = svm.margin, mapping = aes(x = x, y = prob, linetype = method, color = cat)) +
facet_wrap(~var, nrow = 2) + theme_bw() +
labs(x = 'Input Value', y = 'Conditional Probability', linetype = 'SVM Kernel', color = '',
subtitle = 'Utopia Distance + Priority') +
scale_color_manual(labels = c('tru' = 'Expected Marginal', 'start' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'start' = 'skyblue2',
'adapt' = 'red', 'rando' = 'green'),
breaks = c('tru', 'start', 'adapt', 'rando')) +
guides(color = guide_legend(override.aes = list(linetype = c(3, 1, 1, 1)))) +
scale_linetype_discrete(labels = c('SVM-pol' = 'Polynomial', 'rad' = 'SVM-Radial')) +
scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(limits = c(0, 1))
ggplot() +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob - psd, color = cat), linetype = 3) +
geom_path(data = filter(Infer.plt, cat == 'tru'), mapping = aes(x = x, y = prob + psd, color = cat), linetype = 3) +
geom_path(data = filter(Infer.plt, cat != 'tru'), mapping = aes(x = x, y = prob, color = cat)) +
facet_wrap(var~., nrow = 2) +
scale_color_manual(labels = c('tru' = 'Expected Marginal', 'start' = 'Starting Dataset',
'adapt' = '+ Adaptive Sampling', 'rando' = '+ Random Sampling'),
values = c('tru' = 'black', 'start' = 'skyblue2',
'adapt' = 'red', 'rando' = 'green'),
breaks = c('tru', 'start', 'adapt', 'rando')) +
guides(color = guide_legend(override.aes = list(linetype = c(3, 1, 1, 1)))) +
labs(x = '', y = 'Probability of Acceptance', subtitle = 'Cutoff Threshold', color = '') +
scale_y_continuous(expand = c(0, 0.05)) + scale_x_continuous(expand = c(0, 0)) +
theme_bw()
# Calculate coefficients of determination: easier comparison
coefdet = data.frame(method = c(rep('SVM-pol', 3), rep('SVM-rad', 3), rep('GP', 3)),
dataset = c('0start', '1adapt', '2rando'),
coef = c(cor(x = filter(svm.margin, method == 'SVM-pol', cat == 'start')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-pol', cat == 'adapt')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-pol', cat == 'rando')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-rad', cat == 'start')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-rad', cat == 'adapt')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(svm.margin, method == 'SVM-rad', cat == 'rando')$prob,
y = filter(Infer.plt, cat == 'tru')$prob, method = 'pearson')^2,
cor(x = filter(Infer.plt, cat == 'start')$prob,
y = c(filter(Infer.plt, cat == 'tru', var == 'x1')$prob, filter(Infer.plt, cat == 'tru', var == 'x2')$prob),
method = 'pearson')^2,
cor(x = filter(Infer.plt, cat == 'adapt')$prob,
y = c(filter(Infer.plt, cat == 'tru', var == 'x1')$prob, filter(Infer.plt, cat == 'tru', var == 'x2')$prob),
method = 'pearson')^2,
cor(x = filter(Infer.plt, cat == 'rando')$prob,
y = c(filter(Infer.plt, cat == 'tru', var == 'x1')$prob, filter(Infer.plt, cat == 'tru', var == 'x2')$prob),
method = 'pearson')^2))
coefdet
write.csv(coefdet, 'Marginals_CoefDet_radan.csv')
ggplot(coefdet) +
geom_col(mapping = aes(x = dataset, fill = dataset, y = coef))+
facet_grid(~method) +
labs(x = '', y = '1-Variable Marginal Coefficient of Determination', subtitle = 'Utopia Distance + Priority') +
scale_x_discrete(breaks = c()) +
scale_fill_manual(labels = c('0start' = 'Starting Dataset', '1adapt' = '+ Adaptive Sampling',
'2rando' = '+ Random Sampling'), name = '',
values = c('0start' = 'skyblue2', '1adapt' = 'red', '2rando' = 'green')) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05)))
import.marginal = function(svm.margin, typ){
# Set up output
import.svm.margin = data.frame()
# Loop across methods
for(met in unique(svm.margin$method)){
sub = filter(svm.margin, cat == 'adapt', var == 'x1', method == met)
import.svm.margin = rbind(import.svm.margin, data.frame(
import = diff(range(sub$prob))/sum(sub$psd),
var = 'x1',
method = paste('Marginal', substring(met, 4), sep = ''),
sd = sqrt((max(sub$prob)*(1-max(sub$prob)) + min(sub$prob)*(1-min(sub$prob)) *
diff(range(sub$prob)) + var(sub$psd))/nrow(sub))
))
sub = filter(svm.margin, cat == 'adapt', var == 'x2', method == met)
import.svm.margin = rbind(import.svm.margin, data.frame(
import = diff(range(sub$prob))/sum(sub$psd),
var = 'x2',
method = paste('Marginal', substring(met, 4), sep = ''),
sd = sqrt((max(sub$prob)*(1-max(sub$prob)) + min(sub$prob)*(1-min(sub$prob)) *
diff(range(sub$prob)) + var(sub$psd))/nrow(sub))
))
}
# Type, relative importance
import.svm.margin$typ = typ
import.svm.margin$r.import = c(import.svm.margin$import[1:2]/max(import.svm.margin$import[1:2]),
import.svm.margin$import[3:4]/max(import.svm.margin$import[3:4]))
return(import.svm.margin)
}
svm.import.margin = rbind(import.marginal(svm.margin = read.csv('Margin_SVM_dist.csv'), typ = 'Pareto Distance'),
import.marginal(svm.margin = read.csv('Margin_SVM_cutof.csv'), typ = 'Threshold Cutoff'),
import.marginal(svm.margin = read.csv('Margin_SVM_radan.csv'), typ = 'Utopia Distance'))
import.rank = read.csv('../Ex_Quartic/Importance.csv')
import.rank = import.rank[, !names(import.rank) %in% c('X')]
import.rank$method = unlist(lapply(import.rank$method, function(x) paste(x, '-GP', sep = '')))
import.rank
# import.rank$method[import.rank$method == 'Shapley'] = 'Shapley-GP'
SVM.shap.delta = SVM.shap(mod = mod.adapt.delta)
SVM.shap.delta$typ = 'Pareto Distance'
SVM.shap.cutof = SVM.shap(mod = mod.adapt.cutof)
SVM.shap.cutof$typ = 'Threshold Cutoff'
SVM.shap.radan = SVM.shap(mod = mod.adapt.radan)
SVM.shap.radan$typ = 'Utopia Distance'
ggplot(rbind(import.rank, SVM.shap.delta, SVM.shap.cutof, SVM.shap.radan, svm.import.margin)) +
geom_col(mapping = aes(x = var, y = abs(r.import), fill = var)) +
facet_grid(method~typ, scales = 'free_y') +
labs(x = '', y = 'Relative Importance') +
guides(fill = FALSE) + scale_x_discrete(labels = c('x1' = expression('x'[1]), 'x2' = expression('x'[2])))
ggplot(rbind(import.rank, SVM.shap.delta, SVM.shap.cutof, SVM.shap.radan, svm.import.margin)) +
geom_col(mapping = aes(x = var, y = import, fill = var)) +
facet_grid(method~typ, scales = 'free_y') +
labs(x = '', y = 'Relative Importance') +
guides(fill = FALSE) + scale_x_discrete(labels = c('x1' = expression('x'[1]), 'x2' = expression('x'[2])))
# Rows: Selection method = Pareto distance, Objective Cutoff, Utopia Distance + Priority
# Columns: method = Shapley with SVM, Shapley with GP, new method
# Split up by model method
import.all = rbind(import.rank, SVM.shap.delta, SVM.shap.cutof, SVM.shap.radan, svm.import.margin)
# Classification split
split = unlist(lapply(import.all$method, strsplit, '-'))
import.all$method = split[c(TRUE, FALSE)]
import.all$model = split[c(FALSE, TRUE)]
import.all
ggplot(import.all) +
geom_col(mapping = aes(x = model, y = import, fill = var), position = 'dodge') +
geom_errorbar(mapping = aes(x = model, ymax = import + sd, ymin = import - sd, group = var),
position = 'dodge', width = 0.9) +
facet_grid(method~typ, scales = 'free_y') +
labs(x = 'Model', y = 'Relative Importance (Unitless)') +
scale_x_discrete(labels = c('GP', 'pol' = 'SVM-pol', 'rad' = 'SVM-rad')) +
scale_fill_discrete(name = '', labels = c('x1' = expression('x'[1]), 'x2' = expression('x'[2])))
# Relative importance (normalized): uncertainty is divided by the same value
r.max = import.all$import
for(i in seq(from = 1, to = length(r.max)-1, by = 2)){
r.max[c(i,i+1)] = max(abs(r.max[c(i,i+1)]))
}
import.all$r.import = abs(import.all$import)/r.max
import.all$r.sd = import.all$sd / r.max
ggplot(import.all) +
geom_col(mapping = aes(x = model, y = r.import, fill = var), position = 'dodge') +
geom_errorbar(mapping = aes(x = model, ymax = r.import + r.sd, ymin = r.import - r.sd, group = var),
position = 'dodge', width = 0.9) +
facet_grid(method~typ, scales = 'free_y') +
labs(x = 'Model', y = 'Relative Importance (Unitless)') +
scale_x_discrete(labels = c('GP', 'pol' = 'SVM-pol', 'rad' = 'SVM-rad')) +
scale_fill_discrete(name = '', labels = c('x1' = expression('x'[1]), 'x2' = expression('x'[2])))
rm(r.max)
write.csv(import.all, 'ImportRank_Method.csv')
Since the GM models are unsupervised, there cannot be control for what the selection criteria are. A superficial analysis will be conducted to check what the ML algorithm is converging to, but the results will be interpreted solely in terms of descriptive statistics.
For the purposes of analysis, three models will be tested based on the input to the ML algorithm: * (x1, x2, f1, f2): the full set of inputs and outputs - since it has the most data, this will likely be the most robust * (x1, x2): negative control of just the inputs - it should group the points based on the sampling densities * (f1, f2): outputs only - this should have the best distinction between good and bad performing groups
Up to 12 categories will be allowed, and any type of Gaussian models will be accepted. This will lead to longer fitting time, but should help with accuracy.
GPar.all = read.csv(file = '../Ex_Quartic/GPar_all_start.csv')
# names(GPar.all)
max.cat = 12
cluster.all = Mclust(data = GPar.all[,c('x1', 'x2', 'f1', 'f2')], G = 2:max.cat)
fitting ...
|
| | 0%
|
|= | 1%
|
|== | 2%
|
|=== | 3%
|
|==== | 3%
|
|==== | 4%
|
|===== | 5%
|
|====== | 5%
|
|====== | 6%
|
|======= | 6%
|
|======== | 7%
|
|======== | 8%
|
|========= | 8%
|
|========== | 9%
|
|=========== | 10%
|
|============ | 11%
|
|============= | 12%
|
|============== | 13%
|
|=============== | 14%
|
|================ | 15%
|
|================= | 15%
|
|================== | 16%
|
|================== | 17%
|
|=================== | 17%
|
|==================== | 18%
|
|==================== | 19%
|
|===================== | 19%
|
|====================== | 20%
|
|======================= | 21%
|
|======================== | 22%
|
|========================= | 23%
|
|========================== | 24%
|
|=========================== | 25%
|
|============================ | 26%
|
|============================= | 26%
|
|============================== | 27%
|
|============================== | 28%
|
|=============================== | 28%
|
|================================ | 29%
|
|================================ | 30%
|
|================================= | 30%
|
|================================== | 31%
|
|================================== | 32%
|
|=================================== | 32%
|
|==================================== | 33%
|
|===================================== | 34%
|
|====================================== | 35%
|
|======================================= | 35%
|
|======================================= | 36%
|
|======================================== | 37%
|
|========================================= | 37%
|
|========================================= | 38%
|
|========================================== | 39%
|
|=========================================== | 39%
|
|============================================ | 40%
|
|============================================ | 41%
|
|============================================= | 41%
|
|============================================== | 42%
|
|============================================== | 43%
|
|=============================================== | 43%
|
|================================================ | 44%
|
|================================================= | 45%
|
|================================================== | 46%
|
|=================================================== | 46%
|
|=================================================== | 47%
|
|==================================================== | 48%
|
|===================================================== | 48%
|
|===================================================== | 49%
|
|====================================================== | 50%
|
|======================================================= | 50%
|
|======================================================== | 51%
|
|======================================================== | 52%
|
|========================================================= | 52%
|
|========================================================== | 53%
|
|========================================================== | 54%
|
|=========================================================== | 54%
|
|============================================================ | 55%
|
|============================================================= | 56%
|
|============================================================== | 57%
|
|=============================================================== | 57%
|
|=============================================================== | 58%
|
|================================================================ | 59%
|
|================================================================= | 59%
|
|================================================================= | 60%
|
|================================================================== | 61%
|
|=================================================================== | 61%
|
|==================================================================== | 62%
|
|==================================================================== | 63%
|
|===================================================================== | 63%
|
|====================================================================== | 64%
|
|====================================================================== | 65%
|
|======================================================================= | 65%
|
|======================================================================== | 66%
|
|========================================================================= | 67%
|
|========================================================================== | 68%
|
|=========================================================================== | 68%
|
|=========================================================================== | 69%
|
|============================================================================ | 70%
|
|============================================================================= | 70%
|
|============================================================================= | 71%
|
|============================================================================== | 72%
|
|=============================================================================== | 72%
|
|=============================================================================== | 73%
|
|================================================================================ | 74%
|
|================================================================================= | 74%
|
|================================================================================== | 75%
|
|=================================================================================== | 76%
|
|==================================================================================== | 77%
|
|===================================================================================== | 78%
|
|====================================================================================== | 79%
|
|======================================================================================= | 80%
|
|======================================================================================== | 81%
|
|========================================================================================= | 81%
|
|========================================================================================= | 82%
|
|========================================================================================== | 83%
|
|=========================================================================================== | 83%
|
|=========================================================================================== | 84%
|
|============================================================================================ | 85%
|
|============================================================================================= | 85%
|
|============================================================================================== | 86%
|
|=============================================================================================== | 87%
|
|================================================================================================ | 88%
|
|================================================================================================= | 89%
|
|================================================================================================== | 90%
|
|=================================================================================================== | 91%
|
|==================================================================================================== | 92%
|
|===================================================================================================== | 92%
|
|===================================================================================================== | 93%
|
|====================================================================================================== | 94%
|
|======================================================================================================= | 94%
|
|======================================================================================================= | 95%
|
|======================================================================================================== | 95%
|
|========================================================================================================= | 96%
|
|========================================================================================================= | 97%
|
|========================================================================================================== | 97%
|
|=========================================================================================================== | 98%
|
|============================================================================================================ | 99%
|
|=============================================================================================================| 100%
cluster.in = Mclust(data = GPar.all[,c('x1', 'x2')], G = 2:max.cat)
fitting ...
|
| | 0%
|
|= | 1%
|
|== | 2%
|
|=== | 3%
|
|==== | 3%
|
|==== | 4%
|
|===== | 5%
|
|====== | 5%
|
|====== | 6%
|
|======= | 6%
|
|======== | 7%
|
|======== | 8%
|
|========= | 8%
|
|========== | 9%
|
|=========== | 10%
|
|============ | 11%
|
|============= | 12%
|
|============== | 13%
|
|=============== | 14%
|
|================ | 15%
|
|================= | 15%
|
|================== | 16%
|
|================== | 17%
|
|=================== | 17%
|
|==================== | 18%
|
|==================== | 19%
|
|===================== | 19%
|
|====================== | 20%
|
|======================= | 21%
|
|======================== | 22%
|
|========================= | 23%
|
|========================== | 24%
|
|=========================== | 25%
|
|============================ | 26%
|
|============================= | 26%
|
|============================== | 27%
|
|============================== | 28%
|
|=============================== | 28%
|
|================================ | 29%
|
|================================ | 30%
|
|================================= | 30%
|
|================================== | 31%
|
|================================== | 32%
|
|=================================== | 32%
|
|==================================== | 33%
|
|===================================== | 34%
|
|====================================== | 35%
|
|======================================= | 35%
|
|======================================= | 36%
|
|======================================== | 37%
|
|========================================= | 37%
|
|========================================= | 38%
|
|========================================== | 39%
|
|=========================================== | 39%
|
|============================================ | 40%
|
|============================================ | 41%
|
|============================================= | 41%
|
|============================================== | 42%
|
|============================================== | 43%
|
|=============================================== | 43%
|
|================================================ | 44%
|
|================================================= | 45%
|
|================================================== | 46%
|
|=================================================== | 46%
|
|=================================================== | 47%
|
|==================================================== | 48%
|
|===================================================== | 48%
|
|===================================================== | 49%
|
|====================================================== | 50%
|
|======================================================= | 50%
|
|======================================================== | 51%
|
|======================================================== | 52%
|
|========================================================= | 52%
|
|========================================================== | 53%
|
|========================================================== | 54%
|
|=========================================================== | 54%
|
|============================================================ | 55%
|
|============================================================= | 56%
|
|============================================================== | 57%
|
|=============================================================== | 57%
|
|=============================================================== | 58%
|
|================================================================ | 59%
|
|================================================================= | 59%
|
|================================================================= | 60%
|
|================================================================== | 61%
|
|=================================================================== | 61%
|
|==================================================================== | 62%
|
|==================================================================== | 63%
|
|===================================================================== | 63%
|
|====================================================================== | 64%
|
|====================================================================== | 65%
|
|======================================================================= | 65%
|
|======================================================================== | 66%
|
|========================================================================= | 67%
|
|========================================================================== | 68%
|
|=========================================================================== | 68%
|
|=========================================================================== | 69%
|
|============================================================================ | 70%
|
|============================================================================= | 70%
|
|============================================================================= | 71%
|
|============================================================================== | 72%
|
|=============================================================================== | 72%
|
|=============================================================================== | 73%
|
|================================================================================ | 74%
|
|================================================================================= | 74%
|
|================================================================================== | 75%
|
|=================================================================================== | 76%
|
|==================================================================================== | 77%
|
|===================================================================================== | 78%
|
|====================================================================================== | 79%
|
|======================================================================================= | 80%
|
|======================================================================================== | 81%
|
|========================================================================================= | 81%
|
|========================================================================================= | 82%
|
|========================================================================================== | 83%
|
|=========================================================================================== | 83%
|
|=========================================================================================== | 84%
|
|============================================================================================ | 85%
|
|============================================================================================= | 85%
|
|============================================================================================== | 86%
|
|=============================================================================================== | 87%
|
|================================================================================================ | 88%
|
|================================================================================================= | 89%
|
|================================================================================================== | 90%
|
|=================================================================================================== | 91%
|
|==================================================================================================== | 92%
|
|===================================================================================================== | 92%
|
|===================================================================================================== | 93%
|
|====================================================================================================== | 94%
|
|======================================================================================================= | 94%
|
|======================================================================================================= | 95%
|
|======================================================================================================== | 95%
|
|========================================================================================================= | 96%
|
|========================================================================================================= | 97%
|
|========================================================================================================== | 97%
|
|=========================================================================================================== | 98%
|
|============================================================================================================ | 99%
|
|=============================================================================================================| 100%
cluster.out = Mclust(data = GPar.all[,c('f1', 'f2')], G = 2:max.cat)
fitting ...
|
| | 0%
|
|= | 1%
|
|== | 2%
|
|=== | 3%
|
|==== | 3%
|
|==== | 4%
|
|===== | 5%
|
|====== | 5%
|
|====== | 6%
|
|======= | 6%
|
|======== | 7%
|
|======== | 8%
|
|========= | 8%
|
|========== | 9%
|
|=========== | 10%
|
|============ | 11%
|
|============= | 12%
|
|============== | 13%
|
|=============== | 14%
|
|================ | 15%
|
|================= | 15%
|
|================== | 16%
|
|================== | 17%
|
|=================== | 17%
|
|==================== | 18%
|
|==================== | 19%
|
|===================== | 19%
|
|====================== | 20%
|
|======================= | 21%
|
|======================== | 22%
|
|========================= | 23%
|
|========================== | 24%
|
|=========================== | 25%
|
|============================ | 26%
|
|============================= | 26%
|
|============================== | 27%
|
|============================== | 28%
|
|=============================== | 28%
|
|================================ | 29%
|
|================================ | 30%
|
|================================= | 30%
|
|================================== | 31%
|
|================================== | 32%
|
|=================================== | 32%
|
|==================================== | 33%
|
|===================================== | 34%
|
|====================================== | 35%
|
|======================================= | 35%
|
|======================================= | 36%
|
|======================================== | 37%
|
|========================================= | 37%
|
|========================================= | 38%
|
|========================================== | 39%
|
|=========================================== | 39%
|
|============================================ | 40%
|
|============================================ | 41%
|
|============================================= | 41%
|
|============================================== | 42%
|
|============================================== | 43%
|
|=============================================== | 43%
|
|================================================ | 44%
|
|================================================= | 45%
|
|================================================== | 46%
|
|=================================================== | 46%
|
|=================================================== | 47%
|
|==================================================== | 48%
|
|===================================================== | 48%
|
|===================================================== | 49%
|
|====================================================== | 50%
|
|======================================================= | 50%
|
|======================================================== | 51%
|
|======================================================== | 52%
|
|========================================================= | 52%
|
|========================================================== | 53%
|
|========================================================== | 54%
|
|=========================================================== | 54%
|
|============================================================ | 55%
|
|============================================================= | 56%
|
|============================================================== | 57%
|
|=============================================================== | 57%
|
|=============================================================== | 58%
|
|================================================================ | 59%
|
|================================================================= | 59%
|
|================================================================= | 60%
|
|================================================================== | 61%
|
|=================================================================== | 61%
|
|==================================================================== | 62%
|
|==================================================================== | 63%
|
|===================================================================== | 63%
|
|====================================================================== | 64%
|
|====================================================================== | 65%
|
|======================================================================= | 65%
|
|======================================================================== | 66%
|
|========================================================================= | 67%
|
|========================================================================== | 68%
|
|=========================================================================== | 68%
|
|=========================================================================== | 69%
|
|============================================================================ | 70%
|
|============================================================================= | 70%
|
|============================================================================= | 71%
|
|============================================================================== | 72%
|
|=============================================================================== | 72%
|
|=============================================================================== | 73%
|
|================================================================================ | 74%
|
|================================================================================= | 74%
|
|================================================================================== | 75%
|
|=================================================================================== | 76%
|
|==================================================================================== | 77%
|
|===================================================================================== | 78%
|
|====================================================================================== | 79%
|
|======================================================================================= | 80%
|
|======================================================================================== | 81%
|
|========================================================================================= | 81%
|
|========================================================================================= | 82%
|
|========================================================================================== | 83%
|
|=========================================================================================== | 83%
|
|=========================================================================================== | 84%
|
|============================================================================================ | 85%
|
|============================================================================================= | 85%
|
|============================================================================================== | 86%
|
|=============================================================================================== | 87%
|
|================================================================================================ | 88%
|
|================================================================================================= | 89%
|
|================================================================================================== | 90%
|
|=================================================================================================== | 91%
|
|==================================================================================================== | 92%
|
|===================================================================================================== | 92%
|
|===================================================================================================== | 93%
|
|====================================================================================================== | 94%
|
|======================================================================================================= | 94%
|
|======================================================================================================= | 95%
|
|======================================================================================================== | 95%
|
|========================================================================================================= | 96%
|
|========================================================================================================= | 97%
|
|========================================================================================================== | 97%
|
|=========================================================================================================== | 98%
|
|============================================================================================================ | 99%
|
|=============================================================================================================| 100%
# summary(cluster, parameters = TRUE)
GPar.all$typ.all = cluster.all$classification
GPar.all$typ.in = cluster.in$classification
GPar.all$typ.out = cluster.out$classification
ggplot() +
# Boundaries: +/- some separation from 0.5
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = dist, color = 'delta'), breaks = c(1)) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = cutof, color = 'cutof'), breaks = c(0.5)) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = rad, color = 'rad'), breaks = c(1)) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = ang, color = 'ang'), breaks = c(50)) +
# Pareto frontier
geom_smooth(data = GPar.front, mapping = aes(x = x1, y = x2, color = 'Pareto'),
level = 0.95, formula = (y~x), method = 'loess') +
# Clustering
geom_point(data = GPar.all, mapping = aes(x = x1, y = x2, fill = as.factor(typ.all)), size = 2.5, shape = 21) +
labs(x = expression('x'[1]), y = expression('x'[2]),
color = 'Acceptance Criteria', fill = 'Group',
subtitle = '(x1, x2, f1, f2)') +
scale_color_manual(values = c('delta' = '#1b9e77', 'cutof' = '#d95f02',
'rad' = '#7570b3', 'ang' = '#e7298a',
'Pareto' = 'black'),
labels = c('delta' = expression(delta*' < 1'),
'cutof' = expression('F'[1]^'*'*'< 1, F'[2]^'*'*'< 1'),
'rad' = expression('r < 1'),
'ang' = expression('F'[1]*'and F'[2]*' Balance'),
'Pareto' = expression('Pareto Front')),
breaks = c('delta', 'cutof', 'rad', 'ang', 'Pareto')) +
theme_classic() + #theme(legend.position = c(0.85, 0.75)) +
scale_x_continuous(expand = c(0, 0), limits = c(0, 5)) + scale_y_continuous(expand = c(0, 0), limits = c(0, 5)) +
guides(colour = guide_legend(override.aes = list(fill = alpha('white', 1))))
ggplot() +
# Boundaries: +/- some separation from 0.5
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = dist, color = 'delta'), breaks = c(1)) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = cutof, color = 'cutof'), breaks = c(0.5)) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = rad, color = 'rad'), breaks = c(1)) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = ang, color = 'ang'), breaks = c(50)) +
# Pareto frontier
geom_smooth(data = GPar.front, mapping = aes(x = x1, y = x2, color = 'Pareto'),
level = 0.95, formula = (y~x), method = 'loess') +
# Clustering
geom_point(data = GPar.all, mapping = aes(x = x1, y = x2, fill = as.factor(typ.in)), size = 2.5, shape = 21) +
labs(x = expression('x'[1]), y = expression('x'[2]),
color = 'Acceptance Criteria', fill = 'Group',
subtitle = '(x1, x2)') +
scale_color_manual(values = c('delta' = '#1b9e77', 'cutof' = '#d95f02',
'rad' = '#7570b3', 'ang' = '#e7298a',
'Pareto' = 'black'),
labels = c('delta' = expression(delta*' < 1'),
'cutof' = expression('F'[1]^'*'*'< 1, F'[2]^'*'*'< 1'),
'rad' = expression('r < 1'),
'ang' = expression('F'[1]*'and F'[2]*' Balance'),
'Pareto' = expression('Pareto Front')),
breaks = c('delta', 'cutof', 'rad', 'ang', 'Pareto')) +
theme_classic() + #theme(legend.position = c(0.85, 0.75)) +
scale_x_continuous(expand = c(0, 0), limits = c(0, 5)) + scale_y_continuous(expand = c(0, 0), limits = c(0, 5)) +
guides(colour = guide_legend(override.aes = list(fill = alpha('white', 1))))
ggplot() +
# Boundaries: +/- some separation from 0.5
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = dist, color = 'delta'), breaks = c(1)) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = cutof, color = 'cutof'), breaks = c(0.5)) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = rad, color = 'rad'), breaks = c(1)) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = ang, color = 'ang'), breaks = c(50)) +
# Pareto frontier
geom_smooth(data = GPar.front, mapping = aes(x = x1, y = x2, color = 'Pareto'),
level = 0.95, formula = (y~x), method = 'loess') +
# Clustering
geom_point(data = GPar.all, mapping = aes(x = x1, y = x2, fill = as.factor(typ.out)), size = 2.5, shape = 21) +
labs(x = expression('x'[1]), y = expression('x'[2]),
color = 'Acceptance Criteria', fill = 'Group',
subtitle = '(f1, f2)') +
scale_color_manual(values = c('delta' = '#1b9e77', 'cutof' = '#d95f02',
'rad' = '#7570b3', 'ang' = '#e7298a',
'Pareto' = 'black'),
labels = c('delta' = expression(delta*' < 1'),
'cutof' = expression('F'[1]^'*'*'< 1, F'[2]^'*'*'< 1'),
'rad' = expression('r < 1'),
'ang' = expression('F'[1]*'and F'[2]*' Balance'),
'Pareto' = expression('Pareto Front')),
breaks = c('delta', 'cutof', 'rad', 'ang', 'Pareto')) +
theme_classic() + #theme(legend.position = c(0.85, 0.75)) +
scale_x_continuous(expand = c(0, 0), limits = c(0, 5)) + scale_y_continuous(expand = c(0, 0), limits = c(0, 5)) +
guides(colour = guide_legend(override.aes = list(fill = alpha('white', 1))))
The complete input (x1, x2, f1, f2) overcomplicates the system, providing 8 groups. One of these groups is clearly the Pareto frontier, but it does not include any points that are of similar performance. There are points to either side of it suggesting similar performance, but how far they extend is difficult to interpret.
The negative control (x1, x2) gives the expected result of largely grouping based on sampling density. This means the highly sampled Pareto frontier and local minimum are their own groups, and the rest of the space is divided based around the boundary between the highly sampled regions.
The output-only model gives the most useful results, as it groups the Pareto frontier inside of another high-performing group, which also includes the local optimum and the space spanning to it. It roughly lines up with the Pareto distance or threshold criteria; it does not appear to prioritize the angle or utopia distance. Showing only this model at finer resolution as it is the only relevant performing model
res = predict(object = cluster.out, newdata = fine.grid[c('f1', 'f2')])
fine.grid$cluster.out = res$classification
ggplot() +
# Cluster results
geom_point(data = fine.grid, mapping = aes(x = x1, y = x2, color = as.factor(cluster.out))) +
# Boundaries: +/- some separation from 0.5
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = dist), color = 'black', breaks = c(1)) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = cutof), color = 'red', breaks = c(0.5)) +
labs(x = expression('x'[1]), y = expression('x'[2]),
color = 'Group', fill = 'Group',
subtitle = '(f1, f2)') +
theme_classic() + #theme(legend.position = c(0.85, 0.75)) +
scale_x_continuous(expand = c(0, 0), limits = c(0, 5)) + scale_y_continuous(expand = c(0, 0), limits = c(0, 5)) +
guides(colour = guide_legend(override.aes = list(fill = alpha('white', 1))))
res = predict(object = cluster.all, newdata = fine.grid[c('x1', 'x2', 'f1', 'f2')])
fine.grid$cluster.all = res$classification
ggplot() +
# Cluster results
geom_point(data = fine.grid, mapping = aes(x = x1, y = x2, color = as.factor(cluster.all)), size = 4) +
# Boundaries: +/- some separation from 0.5
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = dist), color = 'black', breaks = c(1)) +
geom_contour(data = fine.grid, mapping = aes(x = x1, y = x2, z = cutof), color = 'red', breaks = c(0.5)) +
labs(x = expression('x'[1]), y = expression('x'[2]),
color = 'Group', fill = 'Group',
subtitle = '(x1, x2, f1, f2)') +
theme_classic() + #theme(legend.position = c(0.85, 0.75)) +
scale_x_continuous(expand = c(0, 0), limits = c(0, 5)) + scale_y_continuous(expand = c(0, 0), limits = c(0, 5)) +
guides(colour = guide_legend(override.aes = list(fill = alpha('white', 1))))
ggplot() +
# Cluster results
geom_point(data = filter(fine.grid, f2 < 50, f1 < 200), mapping = aes(x = f1, y = f2, color = as.factor(cluster.all)), size = 2) +
labs(x = expression('f'[1]), y = expression('f'[2]),
color = 'Group', fill = 'Group',
subtitle = '(x1, x2, f1, f2)') +
theme_classic() +
guides(colour = guide_legend(override.aes = list(fill = alpha('white', 1))))
The (f1, f2) model appears to give a cluster that is somewhere between a threshold cutoff and the Pareto distance cutoff. In contrast, the behavior of the (x1, x2, f1, f2) model provides many more groups, including splitting the Pareto front into about three different categories with no obvious analog to why the boundaries are where they are. It appears to roughly fit the same boundaries of Pareto distance or threshold cutoffs, but not very well. A rough approximation is that groups (5, 6) are the Pareto front, groups (2, 3, 7) make up the region close to the Pareto front, and (1, 4, 8) are the region far from the front.
The model for (f1, f2) is going to give the same conditional probabilities as the Pareto distance or threshold acceptance criteria based on this similarity in the shape of the boundary. The interesting result to interpret is the (x1, x2, f1, f2), particularly when marginalizing to (x1, x2). Assuming that calculating (f1, f2) is expensive, the best approximation is that found from the GP models used to find the Pareto front itself. These will give (f1, f2) as a bivariate Gaussian distribution, which can be sampled from to estimate the likelihood that it falls into the Pareto front, the region close to it, or the region far from it. This is achievable with a sequential Monte Carlo: given x1, sample x2 from the range, find the distribution of (f1, f2), sample (f1, f2), and solve the classification into the three groups.
f1.mod = fill.sample.mod(GPar.data = GPar.all, input.name = c('x1', 'x2'), output.name = 'f1')
f2.mod = fill.sample.mod(GPar.data = GPar.all, input.name = c('x1', 'x2'), output.name = 'f2')
x.rng = seq(from = 0, to = 5, length.out = 50)
nsamp.x = 100; nsamp.f = 100
gm.margin = data.frame()
for(x in x.rng){
# x1
# Sample x2
inframe = data.frame(x1 = x, x2 = runif(n = nsamp.x, min = 0, max = 5))
classes = c()
for(n in 1:nrow(inframe)){
# Estimate distribution for f1, f2
f1.est = predict(object = f1.mod, newdata = inframe, type = 'UK')
f2.est = predict(object = f2.mod, newdata = inframe, type = 'UK')
test.frame = data.frame(x1 = inframe$x1, x2 = inframe$x2,
f1 = rnorm(n = nsamp.f, mean = f1.est$mean, sd = f1.est$sd),
f2 = rnorm(n = nsamp.f, mean = f2.est$mean, sd = f2.est$sd))
res = predict(object = cluster.all, newdata = test.frame)
classes = c(classes, res$classification)
}
classes = data.frame(group = classes)
gm.margin = rbind(gm.margin, data.frame(x = x, var = 'x1',
p.pare = nrow(filter(classes, group == 5 | group == 6))/(nsamp.f*nsamp.x),
p.near = nrow(filter(classes, group == 2 | group == 3 | group == 7))/(nsamp.f*nsamp.x),
p.dist = nrow(filter(classes, group == 1 | group == 4 | group == 8))/(nsamp.f*nsamp.x)))
# x2
# Sample x2
inframe = data.frame(x2 = x, x1 = runif(n = nsamp.x, min = 0, max = 5))
classes = c()
for(n in 1:nrow(inframe)){
# Estimate distribution for f1, f2
f1.est = predict(object = f1.mod, newdata = inframe, type = 'UK')
f2.est = predict(object = f2.mod, newdata = inframe, type = 'UK')
test.frame = data.frame(x1 = inframe$x1, x2 = inframe$x2,
f1 = rnorm(n = nsamp.f, mean = f1.est$mean, sd = f1.est$sd),
f2 = rnorm(n = nsamp.f, mean = f2.est$mean, sd = f2.est$sd))
res = predict(object = cluster.all, newdata = test.frame)
classes = c(classes, res$classification)
}
classes = data.frame(group = classes)
gm.margin = rbind(gm.margin, data.frame(x = x, var = 'x2',
p.pare = nrow(filter(classes, group == 5 | group == 6))/(nsamp.f*nsamp.x),
p.near = nrow(filter(classes, group == 2 | group == 3 | group == 7))/(nsamp.f*nsamp.x),
p.dist = nrow(filter(classes, group == 1 | group == 4 | group == 8))/(nsamp.f*nsamp.x)))
}
write.csv(gm.margin, 'Margin_GM.csv')
gm.margin = read.csv('Margin_GM.csv')
gm.margin = data.frame(x = rep(gm.margin$x, 3),
var = rep(gm.margin$var, 3),
prob = c(gm.margin$p.pare, gm.margin$p.near, gm.margin$p.dist),
typ = c(rep('Pareto', nrow(gm.margin)),
rep('Near-Pareto', nrow(gm.margin)),
rep('Suboptimal', nrow(gm.margin))))
ggplot(gm.margin) +
geom_path(mapping = aes(x = x, y = prob, color = typ)) +
facet_wrap(~var, nrow = 2) +
scale_color_brewer(palette = 'Dark2') +
labs(x = 'Input Value', y = 'Conditional Probability', color = 'Region')
gm.margin = read.csv('Margin_GM.csv')
gm.margin = data.frame(x = rep(gm.margin$x),
var = rep(gm.margin$var),
prob = gm.margin$p.pare + gm.margin$p.near,
typ = 'Optimal')
ggplot(gm.margin) +
geom_path(mapping = aes(x = x, y = prob, color = typ)) +
facet_wrap(~var, nrow = 2) +
scale_color_brewer(palette = 'Dark2') +
labs(x = 'Input Value', y = 'Conditional Probability', color = 'Region')
The estimated probability of being optimal (Pareto front group or the near-Pareto group) is similar to that calculated through the other methods for the Pareto distance or threshold criteria. However, because of the numerous variables required to generate the model, a larger and more complicated sampling procedure is necessary for accurate estimates.